Remove old autovect-branch by moving to "dead" directory.
[official-gcc.git] / old-autovect-branch / gcc / testsuite / gfortran.dg / g77 / 20010519-1.f
blobe9336f1b6ab68e51483ce24a3c8155b5647e890d
1 c { dg-do compile }
2 CHARMM Element source/dimb/nmdimb.src 1.1
3 C.##IF DIMB
4 SUBROUTINE NMDIMB(X,Y,Z,NAT3,BNBND,BIMAG,LNOMA,AMASS,DDS,DDSCR,
5 1 PARDDV,DDV,DDM,PARDDF,DDF,PARDDE,DDEV,DD1BLK,
6 2 DD1BLL,NADD,LRAISE,DD1CMP,INBCMP,JNBCMP,
7 3 NPAR,ATMPAR,ATMPAS,BLATOM,PARDIM,NFREG,NFRET,
8 4 PARFRQ,CUTF1,ITMX,TOLDIM,IUNMOD,IUNRMD,
9 5 LBIG,LSCI,ATMPAD,SAVF,NBOND,IB,JB,DDVALM)
10 C-----------------------------------------------------------------------
11 C 01-Jul-1992 David Perahia, Liliane Mouawad
12 C 15-Dec-1994 Herman van Vlijmen
14 C This is the main routine for the mixed-basis diagonalization.
15 C See: L.Mouawad and D.Perahia, Biopolymers (1993), 33, 599,
16 C and: D.Perahia and L.Mouawad, Comput. Chem. (1995), 19, 241.
17 C The method iteratively solves the diagonalization of the
18 C Hessian matrix. To save memory space, it uses a compressed
19 C form of the Hessian, which only contains the nonzero elements.
20 C In the diagonalization process, approximate eigenvectors are
21 C mixed with Cartesian coordinates to form a reduced basis. The
22 C Hessian is then diagonalized in the reduced basis. By iterating
23 C over different sets of Cartesian coordinates the method ultimately
24 C converges to the exact eigenvalues and eigenvectors (up to the
25 C requested accuracy).
26 C If no existing basis set is read, an initial basis will be created
27 C which consists of the low-frequency eigenvectors of diagonal blocks
28 C of the Hessian.
29 C-----------------------------------------------------------------------
30 C-----------------------------------------------------------------------
31 C:::##INCLUDE '~/charmm_fcm/impnon.fcm'
32 C..##IF VAX IRIS HPUX IRIS GNU CSPP OS2 GWS CRAY ALPHA
33 IMPLICIT NONE
34 C..##ENDIF
35 C-----------------------------------------------------------------------
36 C-----------------------------------------------------------------------
37 C:::##INCLUDE '~/charmm_fcm/stream.fcm'
38 LOGICAL LOWER,QLONGL
39 INTEGER MXSTRM,POUTU
40 PARAMETER (MXSTRM=20,POUTU=6)
41 INTEGER NSTRM,ISTRM,JSTRM,OUTU,PRNLEV,WRNLEV,IOLEV
42 COMMON /CASE/ LOWER, QLONGL
43 COMMON /STREAM/ NSTRM,ISTRM,JSTRM(MXSTRM),OUTU,PRNLEV,WRNLEV,IOLEV
44 C..##IF SAVEFCM
45 C..##ENDIF
46 C-----------------------------------------------------------------------
47 C-----------------------------------------------------------------------
48 C:::##INCLUDE '~/charmm_fcm/dimens.fcm'
49 INTEGER LARGE,MEDIUM,SMALL,REDUCE
50 C..##IF QUANTA
51 C..##ELIF T3D
52 C..##ELSE
53 PARAMETER (LARGE=60120, MEDIUM=25140, SMALL=6120)
54 C..##ENDIF
55 PARAMETER (REDUCE=15000)
56 INTEGER SIZE
57 C..##IF XLARGE
58 C..##ELIF XXLARGE
59 C..##ELIF LARGE
60 C..##ELIF MEDIUM
61 PARAMETER (SIZE=MEDIUM)
62 C..##ELIF REDUCE
63 C..##ELIF SMALL
64 C..##ELIF XSMALL
65 C..##ENDIF
66 C..##IF MMFF
67 integer MAXDEFI
68 parameter(MAXDEFI=250)
69 INTEGER NAME0,NAMEQ0,NRES0,KRES0
70 PARAMETER (NAME0=4,NAMEQ0=10,NRES0=4,KRES0=4)
71 integer MaxAtN
72 parameter (MaxAtN=55)
73 INTEGER MAXAUX
74 PARAMETER (MAXAUX = 10)
75 C..##ENDIF
76 INTEGER MAXCSP, MAXHSET
77 C..##IF HMCM
78 PARAMETER (MAXHSET = 200)
79 C..##ELSE
80 C..##ENDIF
81 C..##IF REDUCE
82 C..##ELSE
83 PARAMETER (MAXCSP = 500)
84 C..##ENDIF
85 C..##IF HMCM
86 INTEGER MAXHCM,MAXPCM,MAXRCM
87 C...##IF REDUCE
88 C...##ELSE
89 PARAMETER (MAXHCM=500)
90 PARAMETER (MAXPCM=5000)
91 PARAMETER (MAXRCM=2000)
92 C...##ENDIF
93 C..##ENDIF
94 INTEGER MXCMSZ
95 C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
96 C..##ELSE
97 PARAMETER (MXCMSZ = 5000)
98 C..##ENDIF
99 INTEGER CHRSIZ
100 PARAMETER (CHRSIZ = SIZE)
101 INTEGER MAXATB
102 C..##IF REDUCE
103 C..##ELIF QUANTA
104 C..##ELSE
105 PARAMETER (MAXATB = 200)
106 C..##ENDIF
107 INTEGER MAXVEC
108 C..##IFN VECTOR PARVECT
109 PARAMETER (MAXVEC = 10)
110 C..##ELIF LARGE XLARGE XXLARGE
111 C..##ELIF MEDIUM
112 C..##ELIF SMALL REDUCE
113 C..##ELIF XSMALL
114 C..##ELSE
115 C..##ENDIF
116 INTEGER IATBMX
117 PARAMETER (IATBMX = 8)
118 INTEGER MAXHB
119 C..##IF LARGE XLARGE XXLARGE
120 C..##ELIF MEDIUM
121 PARAMETER (MAXHB = 8000)
122 C..##ELIF SMALL
123 C..##ELIF REDUCE XSMALL
124 C..##ELSE
125 C..##ENDIF
126 INTEGER MAXTRN,MAXSYM
127 C..##IFN NOIMAGES
128 PARAMETER (MAXTRN = 5000)
129 PARAMETER (MAXSYM = 192)
130 C..##ELSE
131 C..##ENDIF
132 C..##IF LONEPAIR (lonepair_max)
133 INTEGER MAXLP,MAXLPH
134 C...##IF REDUCE
135 C...##ELSE
136 PARAMETER (MAXLP = 2000)
137 PARAMETER (MAXLPH = 4000)
138 C...##ENDIF
139 C..##ENDIF (lonepair_max)
140 INTEGER NOEMAX,NOEMX2
141 C..##IF REDUCE
142 C..##ELSE
143 PARAMETER (NOEMAX = 2000)
144 PARAMETER (NOEMX2 = 4000)
145 C..##ENDIF
146 INTEGER MAXATC, MAXCB, MAXCH, MAXCI, MAXCP, MAXCT, MAXITC, MAXNBF
147 C..##IF REDUCE
148 C..##ELIF MMFF CFF
149 PARAMETER (MAXATC = 500, MAXCB = 1500, MAXCH = 3200, MAXCI = 600,
150 & MAXCP = 3000,MAXCT = 15500,MAXITC = 200, MAXNBF=1000)
151 C..##ELIF YAMMP
152 C..##ELIF LARGE
153 C..##ELSE
154 C..##ENDIF
155 INTEGER MAXCN
156 PARAMETER (MAXCN = MAXITC*(MAXITC+1)/2)
157 INTEGER MAXA, MAXAIM, MAXB, MAXT, MAXP
158 INTEGER MAXIMP, MAXNB, MAXPAD, MAXRES
159 INTEGER MAXSEG, MAXGRP
160 C..##IF LARGE XLARGE XXLARGE
161 C..##ELIF MEDIUM
162 PARAMETER (MAXA = SIZE, MAXB = SIZE, MAXT = SIZE,
163 & MAXP = 2*SIZE)
164 PARAMETER (MAXIMP = 9200, MAXNB = 17200, MAXPAD = 8160,
165 & MAXRES = 14000)
166 C...##IF MCSS
167 C...##ELSE
168 PARAMETER (MAXSEG = 1000)
169 C...##ENDIF
170 C..##ELIF SMALL
171 C..##ELIF XSMALL
172 C..##ELIF REDUCE
173 C..##ELSE
174 C..##ENDIF
175 C..##IF NOIMAGES
176 C..##ELSE
177 PARAMETER (MAXAIM = 2*SIZE)
178 PARAMETER (MAXGRP = 2*SIZE/3)
179 C..##ENDIF
180 INTEGER REDMAX,REDMX2
181 C..##IF REDUCE
182 C..##ELSE
183 PARAMETER (REDMAX = 20)
184 PARAMETER (REDMX2 = 80)
185 C..##ENDIF
186 INTEGER MXRTRS, MXRTA, MXRTB, MXRTT, MXRTP, MXRTI, MXRTX,
187 & MXRTHA, MXRTHD, MXRTBL, NICM
188 PARAMETER (MXRTRS = 200, MXRTA = 5000, MXRTB = 5000,
189 & MXRTT = 5000, MXRTP = 5000, MXRTI = 2000,
190 C..##IF YAMMP
191 C..##ELSE
192 & MXRTX = 5000, MXRTHA = 300, MXRTHD = 300,
193 C..##ENDIF
194 & MXRTBL = 5000, NICM = 10)
195 INTEGER NMFTAB, NMCTAB, NMCATM, NSPLIN
196 C..##IF REDUCE
197 C..##ELSE
198 PARAMETER (NMFTAB = 200, NMCTAB = 3, NMCATM = 12000, NSPLIN = 3)
199 C..##ENDIF
200 INTEGER MAXSHK
201 C..##IF XSMALL
202 C..##ELIF REDUCE
203 C..##ELSE
204 PARAMETER (MAXSHK = SIZE*3/4)
205 C..##ENDIF
206 INTEGER SCRMAX
207 C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
208 C..##ELSE
209 PARAMETER (SCRMAX = 5000)
210 C..##ENDIF
211 C..##IF TSM
212 INTEGER MXPIGG
213 C...##IF REDUCE
214 C...##ELSE
215 PARAMETER (MXPIGG=500)
216 C...##ENDIF
217 INTEGER MXCOLO,MXPUMB
218 PARAMETER (MXCOLO=20,MXPUMB=20)
219 C..##ENDIF
220 C..##IF ADUMB
221 INTEGER MAXUMP, MAXEPA, MAXNUM
222 C...##IF REDUCE
223 C...##ELSE
224 PARAMETER (MAXUMP = 10, MAXNUM = 4)
225 C...##ENDIF
226 C..##ENDIF
227 INTEGER MAXING
228 PARAMETER (MAXING=1000)
229 C..##IF MMFF
230 integer MAX_RINGSIZE, MAX_EACH_SIZE
231 parameter (MAX_RINGSIZE = 20, MAX_EACH_SIZE = 1000)
232 integer MAXPATHS
233 parameter (MAXPATHS = 8000)
234 integer MAX_TO_SEARCH
235 parameter (MAX_TO_SEARCH = 6)
236 C..##ENDIF
237 C-----------------------------------------------------------------------
238 C-----------------------------------------------------------------------
239 C:::##INCLUDE '~/charmm_fcm/number.fcm'
240 REAL(KIND=8) ZERO, ONE, TWO, THREE, FOUR, FIVE, SIX,
241 & SEVEN, EIGHT, NINE, TEN, ELEVEN, TWELVE, THIRTN,
242 & FIFTN, NINETN, TWENTY, THIRTY
243 C..##IF SINGLE
244 C..##ELSE
245 PARAMETER (ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0,
246 & THREE = 3.D0, FOUR = 4.D0, FIVE = 5.D0,
247 & SIX = 6.D0, SEVEN = 7.D0, EIGHT = 8.D0,
248 & NINE = 9.D0, TEN = 10.D0, ELEVEN = 11.D0,
249 & TWELVE = 12.D0, THIRTN = 13.D0, FIFTN = 15.D0,
250 & NINETN = 19.D0, TWENTY = 20.D0, THIRTY = 30.D0)
251 C..##ENDIF
252 REAL(KIND=8) FIFTY, SIXTY, SVNTY2, EIGHTY, NINETY, HUNDRD,
253 & ONE2TY, ONE8TY, THRHUN, THR6TY, NINE99, FIFHUN, THOSND,
254 & FTHSND,MEGA
255 C..##IF SINGLE
256 C..##ELSE
257 PARAMETER (FIFTY = 50.D0, SIXTY = 60.D0, SVNTY2 = 72.D0,
258 & EIGHTY = 80.D0, NINETY = 90.D0, HUNDRD = 100.D0,
259 & ONE2TY = 120.D0, ONE8TY = 180.D0, THRHUN = 300.D0,
260 & THR6TY=360.D0, NINE99 = 999.D0, FIFHUN = 1500.D0,
261 & THOSND = 1000.D0,FTHSND = 5000.D0, MEGA = 1.0D6)
262 C..##ENDIF
263 REAL(KIND=8) MINONE, MINTWO, MINSIX
264 PARAMETER (MINONE = -1.D0, MINTWO = -2.D0, MINSIX = -6.D0)
265 REAL(KIND=8) TENM20,TENM14,TENM8,TENM5,PT0001,PT0005,PT001,PT005,
266 & PT01, PT02, PT05, PTONE, PT125, PT25, SIXTH, THIRD,
267 & PTFOUR, PTSIX, HALF, PT75, PT9999, ONEPT5, TWOPT4
268 C..##IF SINGLE
269 C..##ELSE
270 PARAMETER (TENM20 = 1.0D-20, TENM14 = 1.0D-14, TENM8 = 1.0D-8,
271 & TENM5 = 1.0D-5, PT0001 = 1.0D-4, PT0005 = 5.0D-4,
272 & PT001 = 1.0D-3, PT005 = 5.0D-3, PT01 = 0.01D0,
273 & PT02 = 0.02D0, PT05 = 0.05D0, PTONE = 0.1D0,
274 & PT125 = 0.125D0, SIXTH = ONE/SIX,PT25 = 0.25D0,
275 & THIRD = ONE/THREE,PTFOUR = 0.4D0, HALF = 0.5D0,
276 & PTSIX = 0.6D0, PT75 = 0.75D0, PT9999 = 0.9999D0,
277 & ONEPT5 = 1.5D0, TWOPT4 = 2.4D0)
278 C..##ENDIF
279 REAL(KIND=8) ANUM,FMARK
280 REAL(KIND=8) RSMALL,RBIG
281 C..##IF SINGLE
282 C..##ELSE
283 PARAMETER (ANUM=9999.0D0, FMARK=-999.0D0)
284 PARAMETER (RSMALL=1.0D-10,RBIG=1.0D20)
285 C..##ENDIF
286 REAL(KIND=8) RPRECI,RBIGST
287 C..##IF VAX DEC
288 C..##ELIF IBM
289 C..##ELIF CRAY
290 C..##ELIF ALPHA T3D T3E
291 C..##ELSE
292 C...##IF SINGLE
293 C...##ELSE
294 PARAMETER (RPRECI = 2.22045D-16, RBIGST = 4.49423D+307)
295 C...##ENDIF
296 C..##ENDIF
297 C-----------------------------------------------------------------------
298 C-----------------------------------------------------------------------
299 C:::##INCLUDE '~/charmm_fcm/consta.fcm'
300 REAL(KIND=8) PI,RADDEG,DEGRAD,TWOPI
301 PARAMETER(PI=3.141592653589793D0,TWOPI=2.0D0*PI)
302 PARAMETER (RADDEG=180.0D0/PI)
303 PARAMETER (DEGRAD=PI/180.0D0)
304 REAL(KIND=8) COSMAX
305 PARAMETER (COSMAX=0.9999999999D0)
306 REAL(KIND=8) TIMFAC
307 PARAMETER (TIMFAC=4.88882129D-02)
308 REAL(KIND=8) KBOLTZ
309 PARAMETER (KBOLTZ=1.987191D-03)
310 REAL(KIND=8) CCELEC
311 C..##IF AMBER
312 C..##ELIF DISCOVER
313 C..##ELSE
314 PARAMETER (CCELEC=332.0716D0)
315 C..##ENDIF
316 REAL(KIND=8) CNVFRQ
317 PARAMETER (CNVFRQ=2045.5D0/(2.99793D0*6.28319D0))
318 REAL(KIND=8) SPEEDL
319 PARAMETER (SPEEDL=2.99793D-02)
320 REAL(KIND=8) ATMOSP
321 PARAMETER (ATMOSP=1.4584007D-05)
322 REAL(KIND=8) PATMOS
323 PARAMETER (PATMOS = 1.D0 / ATMOSP )
324 REAL(KIND=8) BOHRR
325 PARAMETER (BOHRR = 0.529177249D0 )
326 REAL(KIND=8) TOKCAL
327 PARAMETER (TOKCAL = 627.5095D0 )
328 C..##IF MMFF
329 REAL(KIND=8) MDAKCAL
330 parameter(MDAKCAL=143.9325D0)
331 C..##ENDIF
332 REAL(KIND=8) DEBYEC
333 PARAMETER ( DEBYEC = 2.541766D0 / BOHRR )
334 REAL(KIND=8) ZEROC
335 PARAMETER ( ZEROC = 298.15D0 )
336 C-----------------------------------------------------------------------
337 C-----------------------------------------------------------------------
338 C:::##INCLUDE '~/charmm_fcm/exfunc.fcm'
339 C..##IF ACE
340 C..##ENDIF
341 C..##IF ADUMB
342 C..##ENDIF
343 CHARACTER*4 GTRMA, NEXTA4, CURRA4
344 CHARACTER*6 NEXTA6
345 CHARACTER*8 NEXTA8
346 CHARACTER*20 NEXT20
347 INTEGER ALLCHR, ALLSTK, ALLHP, DECODI, FIND52,
348 * GETATN, GETRES, GETRSN, GETSEG, GTRMI, I4VAL,
349 * ICHAR4, ICMP16, ILOGI4, INDX, INDXA, INDXAF,
350 * INDXRA, INTEG4, IREAL4, IREAL8, LOCDIF,
351 * LUNASS, MATOM, NEXTI, NINDX, NSELCT, NSELCTV, ATMSEL,
352 * PARNUM, PARINS,
353 * SRCHWD, SRCHWS, STRLNG, DSIZE, SSIZE
354 C..##IF ACE
355 * ,GETNNB
356 C..##ENDIF
357 LOGICAL CHKPTR, EQST, EQSTA, EQSTWC, EQWDWC, DOTRIM, CHECQUE,
358 * HYDROG, INITIA, LONE, LTSTEQ, ORDER, ORDER5,
359 * ORDERR, USEDDT, QTOKDEL, QDIGIT, QALPHA
360 REAL(KIND=8) DECODF, DOTVEC, GTRMF, LENVEC, NEXTF, RANDOM, GTRR8,
361 * RANUMB, R8VAL, RETVAL8, SUMVEC
362 C..##IF ADUMB
363 * ,UMFI
364 C..##ENDIF
365 EXTERNAL GTRMA, NEXTA4, CURRA4, NEXTA6, NEXTA8,NEXT20,
366 * ALLCHR, ALLSTK, ALLHP, DECODI, FIND52,
367 * GETATN, GETRES, GETRSN, GETSEG, GTRMI, I4VAL,
368 * ICHAR4, ICMP16, ILOGI4, INDX, INDXA, INDXAF,
369 * INDXRA, INTEG4, IREAL4, IREAL8, LOCDIF,
370 * LUNASS, MATOM, NEXTI, NINDX, NSELCT, NSELCTV, ATMSEL,
371 * PARNUM, PARINS,
372 * SRCHWD, SRCHWS, STRLNG, DSIZE, SSIZE,
373 * CHKPTR, EQST, EQSTA, EQSTWC, EQWDWC, DOTRIM, CHECQUE,
374 * HYDROG, INITIA, LONE, LTSTEQ, ORDER, ORDER5,
375 * ORDERR, USEDDT, QTOKDEL, QDIGIT, QALPHA,
376 * DECODF, DOTVEC, GTRMF, LENVEC, NEXTF, RANDOM, GTRR8,
377 * RANUMB, R8VAL, RETVAL8, SUMVEC
378 C..##IF ADUMB
379 * ,UMFI
380 C..##ENDIF
381 C..##IF ACE
382 * ,GETNNB
383 C..##ENDIF
384 C..##IFN NOIMAGES
385 INTEGER IMATOM
386 EXTERNAL IMATOM
387 C..##ENDIF
388 C..##IF MBOND
389 C..##ENDIF
390 C..##IF MMFF
391 INTEGER LEN_TRIM
392 EXTERNAL LEN_TRIM
393 CHARACTER*4 AtName
394 external AtName
395 CHARACTER*8 ElementName
396 external ElementName
397 CHARACTER*10 QNAME
398 external QNAME
399 integer IATTCH, IBORDR, CONN12, CONN13, CONN14
400 integer LEQUIV, LPATH
401 integer nbndx, nbnd2, nbnd3, NTERMA
402 external IATTCH, IBORDR, CONN12, CONN13, CONN14
403 external LEQUIV, LPATH
404 external nbndx, nbnd2, nbnd3, NTERMA
405 external find_loc
406 REAL(KIND=8) vangle, OOPNGL, TORNGL, ElementMass
407 external vangle, OOPNGL, TORNGL, ElementMass
408 C..##ENDIF
409 C-----------------------------------------------------------------------
410 C-----------------------------------------------------------------------
411 C:::##INCLUDE '~/charmm_fcm/stack.fcm'
412 INTEGER STKSIZ
413 C..##IFN UNICOS
414 C...##IF LARGE XLARGE
415 C...##ELIF MEDIUM REDUCE
416 PARAMETER (STKSIZ=4000000)
417 C...##ELIF SMALL
418 C...##ELIF XSMALL
419 C...##ELIF XXLARGE
420 C...##ELSE
421 C...##ENDIF
422 INTEGER LSTUSD,MAXUSD,STACK
423 COMMON /ISTACK/ LSTUSD,MAXUSD,STACK(STKSIZ)
424 C..##ELSE
425 C..##ENDIF
426 C..##IF SAVEFCM
427 C..##ENDIF
428 C-----------------------------------------------------------------------
429 C-----------------------------------------------------------------------
430 C:::##INCLUDE '~/charmm_fcm/heap.fcm'
431 INTEGER HEAPDM
432 C..##IFN UNICOS (unicos)
433 C...##IF XXLARGE (size)
434 C...##ELIF LARGE XLARGE (size)
435 C...##ELIF MEDIUM (size)
436 C....##IF T3D (t3d2)
437 C....##ELIF TERRA (t3d2)
438 C....##ELIF ALPHA (t3d2)
439 C....##ELIF T3E (t3d2)
440 C....##ELSE (t3d2)
441 PARAMETER (HEAPDM=2048000)
442 C....##ENDIF (t3d2)
443 C...##ELIF SMALL (size)
444 C...##ELIF REDUCE (size)
445 C...##ELIF XSMALL (size)
446 C...##ELSE (size)
447 C...##ENDIF (size)
448 INTEGER FREEHP,HEAPSZ,HEAP
449 COMMON /HEAPST/ FREEHP,HEAPSZ,HEAP(HEAPDM)
450 LOGICAL LHEAP(HEAPDM)
451 EQUIVALENCE (LHEAP,HEAP)
452 C..##ELSE (unicos)
453 C..##ENDIF (unicos)
454 C..##IF SAVEFCM (save)
455 C..##ENDIF (save)
456 C-----------------------------------------------------------------------
457 C-----------------------------------------------------------------------
458 C:::##INCLUDE '~/charmm_fcm/fast.fcm'
459 INTEGER IACNB, NITCC, ICUSED, FASTER, LFAST, LMACH, OLMACH
460 INTEGER ICCOUNT, LOWTP, IGCNB, NITCC2
461 INTEGER ICCNBA, ICCNBB, ICCNBC, ICCNBD, LCCNBA, LCCNBD
462 COMMON /FASTI/ FASTER, LFAST, LMACH, OLMACH, NITCC, NITCC2,
463 & ICUSED(MAXATC), ICCOUNT(MAXATC), LOWTP(MAXATC),
464 & IACNB(MAXAIM), IGCNB(MAXATC),
465 & ICCNBA, ICCNBB, ICCNBC, ICCNBD, LCCNBA, LCCNBD
466 C..##IF SAVEFCM
467 C..##ENDIF
468 C-----------------------------------------------------------------------
469 C-----------------------------------------------------------------------
470 C:::##INCLUDE '~/charmm_fcm/deriv.fcm'
471 REAL(KIND=8) DX,DY,DZ
472 COMMON /DERIVR/ DX(MAXAIM),DY(MAXAIM),DZ(MAXAIM)
473 C..##IF SAVEFCM
474 C..##ENDIF
475 C-----------------------------------------------------------------------
476 C-----------------------------------------------------------------------
477 C:::##INCLUDE '~/charmm_fcm/energy.fcm'
478 INTEGER LENENP, LENENT, LENENV, LENENA
479 PARAMETER (LENENP = 50, LENENT = 70, LENENV = 50,
480 & LENENA = LENENP + LENENT + LENENV )
481 INTEGER TOTE, TOTKE, EPOT, TEMPS, GRMS, BPRESS, PJNK1, PJNK2,
482 & PJNK3, PJNK4, HFCTE, HFCKE, EHFC, EWORK, VOLUME, PRESSE,
483 & PRESSI, VIRI, VIRE, VIRKE, TEPR, PEPR, KEPR, KEPR2,
484 & DROFFA,
485 & XTLTE, XTLKE, XTLPE, XTLTEM, XTLPEP, XTLKEP, XTLKP2,
486 & TOT4, TOTK4, EPOT4, TEM4, MbMom, BodyT, PartT
487 C..##IF ACE
488 & , SELF, SCREEN, COUL ,SOLV, INTER
489 C..##ENDIF
490 C..##IF FLUCQ
491 & ,FQKIN
492 C..##ENDIF
493 PARAMETER (TOTE = 1, TOTKE = 2, EPOT = 3, TEMPS = 4,
494 & GRMS = 5, BPRESS = 6, PJNK1 = 7, PJNK2 = 8,
495 & PJNK3 = 9, PJNK4 = 10, HFCTE = 11, HFCKE = 12,
496 & EHFC = 13, EWORK = 11, VOLUME = 15, PRESSE = 16,
497 & PRESSI = 17, VIRI = 18, VIRE = 19, VIRKE = 20,
498 & TEPR = 21, PEPR = 22, KEPR = 23, KEPR2 = 24,
499 & DROFFA = 26, XTLTE = 27, XTLKE = 28,
500 & XTLPE = 29, XTLTEM = 30, XTLPEP = 31, XTLKEP = 32,
501 & XTLKP2 = 33,
502 & TOT4 = 37, TOTK4 = 38, EPOT4 = 39, TEM4 = 40,
503 & MbMom = 41, BodyT = 42, PartT = 43
504 C..##IF ACE
505 & , SELF = 45, SCREEN = 46, COUL = 47,
506 & SOLV = 48, INTER = 49
507 C..##ENDIF
508 C..##IF FLUCQ
509 & ,FQKIN = 50
510 C..##ENDIF
512 C..##IF ACE
513 C..##ENDIF
514 C..##IF GRID
515 C..##ENDIF
516 C..##IF FLUCQ
517 C..##ENDIF
518 INTEGER BOND, ANGLE, UREYB, DIHE, IMDIHE, VDW, ELEC, HBOND,
519 & USER, CHARM, CDIHE, CINTCR, CQRT, NOE, SBNDRY,
520 & IMVDW, IMELEC, IMHBND, EWKSUM, EWSELF, EXTNDE, RXNFLD,
521 & ST2, IMST2, TSM, QMEL, QMVDW, ASP, EHARM, GEO, MDIP,
522 & PRMS, PANG, SSBP, BK4D, SHEL, RESD, SHAP,
523 & STRB, OOPL, PULL, POLAR, DMC, RGY, EWEXCL, EWQCOR,
524 & EWUTIL, PBELEC, PBNP, PINT, MbDefrm, MbElec, STRSTR,
525 & BNDBND, BNDTW, EBST, MBST, BBT, SST, GBEnr, GSBP
526 C..##IF HMCM
527 & , HMCM
528 C..##ENDIF
529 C..##IF ADUMB
530 & , ADUMB
531 C..##ENDIF
532 & , HYDR
533 C..##IF FLUCQ
534 & , FQPOL
535 C..##ENDIF
536 PARAMETER (BOND = 1, ANGLE = 2, UREYB = 3, DIHE = 4,
537 & IMDIHE = 5, VDW = 6, ELEC = 7, HBOND = 8,
538 & USER = 9, CHARM = 10, CDIHE = 11, CINTCR = 12,
539 & CQRT = 13, NOE = 14, SBNDRY = 15, IMVDW = 16,
540 & IMELEC = 17, IMHBND = 18, EWKSUM = 19, EWSELF = 20,
541 & EXTNDE = 21, RXNFLD = 22, ST2 = 23, IMST2 = 24,
542 & TSM = 25, QMEL = 26, QMVDW = 27, ASP = 28,
543 & EHARM = 29, GEO = 30, MDIP = 31, PINT = 32,
544 & PRMS = 33, PANG = 34, SSBP = 35, BK4D = 36,
545 & SHEL = 37, RESD = 38, SHAP = 39, STRB = 40,
546 & OOPL = 41, PULL = 42, POLAR = 43, DMC = 44,
547 & RGY = 45, EWEXCL = 46, EWQCOR = 47, EWUTIL = 48,
548 & PBELEC = 49, PBNP = 50, MbDefrm= 51, MbElec = 52,
549 & STRSTR = 53, BNDBND = 54, BNDTW = 55, EBST = 56,
550 & MBST = 57, BBT = 58, SST = 59, GBEnr = 60,
551 & GSBP = 65
552 C..##IF HMCM
553 & , HMCM = 61
554 C..##ENDIF
555 C..##IF ADUMB
556 & , ADUMB = 62
557 C..##ENDIF
558 & , HYDR = 63
559 C..##IF FLUCQ
560 & , FQPOL = 65
561 C..##ENDIF
563 INTEGER VEXX, VEXY, VEXZ, VEYX, VEYY, VEYZ, VEZX, VEZY, VEZZ,
564 & VIXX, VIXY, VIXZ, VIYX, VIYY, VIYZ, VIZX, VIZY, VIZZ,
565 & PEXX, PEXY, PEXZ, PEYX, PEYY, PEYZ, PEZX, PEZY, PEZZ,
566 & PIXX, PIXY, PIXZ, PIYX, PIYY, PIYZ, PIZX, PIZY, PIZZ
567 PARAMETER ( VEXX = 1, VEXY = 2, VEXZ = 3, VEYX = 4,
568 & VEYY = 5, VEYZ = 6, VEZX = 7, VEZY = 8,
569 & VEZZ = 9,
570 & VIXX = 10, VIXY = 11, VIXZ = 12, VIYX = 13,
571 & VIYY = 14, VIYZ = 15, VIZX = 16, VIZY = 17,
572 & VIZZ = 18,
573 & PEXX = 19, PEXY = 20, PEXZ = 21, PEYX = 22,
574 & PEYY = 23, PEYZ = 24, PEZX = 25, PEZY = 26,
575 & PEZZ = 27,
576 & PIXX = 28, PIXY = 29, PIXZ = 30, PIYX = 31,
577 & PIYY = 32, PIYZ = 33, PIZX = 34, PIZY = 35,
578 & PIZZ = 36)
579 CHARACTER*4 CEPROP, CETERM, CEPRSS
580 COMMON /ANER/ CEPROP(LENENP), CETERM(LENENT), CEPRSS(LENENV)
581 LOGICAL QEPROP, QETERM, QEPRSS
582 COMMON /QENER/ QEPROP(LENENP), QETERM(LENENT), QEPRSS(LENENV)
583 REAL(KIND=8) EPROP, ETERM, EPRESS
584 COMMON /ENER/ EPROP(LENENP), ETERM(LENENT), EPRESS(LENENV)
585 C..##IF SAVEFCM
586 C..##ENDIF
587 REAL(KIND=8) EPRPA, EPRP2A, EPRPP, EPRP2P,
588 & ETRMA, ETRM2A, ETRMP, ETRM2P,
589 & EPRSA, EPRS2A, EPRSP, EPRS2P
590 COMMON /ENACCM/ EPRPA(LENENP), ETRMA(LENENT), EPRSA(LENENV),
591 & EPRP2A(LENENP),ETRM2A(LENENT),EPRS2A(LENENV),
592 & EPRPP(LENENP), ETRMP(LENENT), EPRSP(LENENV),
593 & EPRP2P(LENENP),ETRM2P(LENENT),EPRS2P(LENENV)
594 C..##IF SAVEFCM
595 C..##ENDIF
596 INTEGER ECALLS, TOT1ST, TOT2ND
597 COMMON /EMISCI/ ECALLS, TOT1ST, TOT2ND
598 REAL(KIND=8) EOLD, FITA, DRIFTA, EAT0A, CORRA, FITP, DRIFTP,
599 & EAT0P, CORRP
600 COMMON /EMISCR/ EOLD, FITA, DRIFTA, EAT0A, CORRA,
601 & FITP, DRIFTP, EAT0P, CORRP
602 C..##IF SAVEFCM
603 C..##ENDIF
604 C..##IF ACE
605 C..##ENDIF
606 C..##IF FLUCQ
607 C..##ENDIF
608 C..##IF ADUMB
609 C..##ENDIF
610 C..##IF GRID
611 C..##ENDIF
612 C..##IF FLUCQ
613 C..##ENDIF
614 C..##IF TSM
615 REAL(KIND=8) TSMTRM(LENENT),TSMTMP(LENENT)
616 COMMON /TSMENG/ TSMTRM,TSMTMP
617 C...##IF SAVEFCM
618 C...##ENDIF
619 C..##ENDIF
620 REAL(KIND=8) EHQBM
621 LOGICAL HQBM
622 COMMON /HQBMVAR/HQBM
623 C..##IF SAVEFCM
624 C..##ENDIF
625 C-----------------------------------------------------------------------
626 C-----------------------------------------------------------------------
627 C:::##INCLUDE '~/charmm_fcm/dimb.fcm'
628 C..##IF DIMB (dimbfcm)
629 INTEGER NPARMX,MNBCMP,LENDSK
630 PARAMETER (NPARMX=1000,MNBCMP=300,LENDSK=200000)
631 INTEGER IJXXCM,IJXYCM,IJXZCM,IJYXCM,IJYYCM
632 INTEGER IJYZCM,IJZXCM,IJZYCM,IJZZCM
633 INTEGER IIXXCM,IIXYCM,IIXZCM,IIYYCM
634 INTEGER IIYZCM,IIZZCM
635 INTEGER JJXXCM,JJXYCM,JJXZCM,JJYYCM
636 INTEGER JJYZCM,JJZZCM
637 PARAMETER (IJXXCM=1,IJXYCM=2,IJXZCM=3,IJYXCM=4,IJYYCM=5)
638 PARAMETER (IJYZCM=6,IJZXCM=7,IJZYCM=8,IJZZCM=9)
639 PARAMETER (IIXXCM=1,IIXYCM=2,IIXZCM=3,IIYYCM=4)
640 PARAMETER (IIYZCM=5,IIZZCM=6)
641 PARAMETER (JJXXCM=1,JJXYCM=2,JJXZCM=3,JJYYCM=4)
642 PARAMETER (JJYZCM=5,JJZZCM=6)
643 INTEGER ITER,IPAR1,IPAR2,NFSAV,PINBCM,PJNBCM,PDD1CM,LENCMP
644 LOGICAL QDISK,QDW,QCMPCT
645 COMMON /DIMBI/ ITER,IPAR1,IPAR2,NFSAV,PINBCM,PJNBCM,LENCMP
646 COMMON /DIMBL/ QDISK,QDW,QCMPCT
647 C...##IF SAVEFCM
648 C...##ENDIF
649 C..##ENDIF (dimbfcm)
650 C-----------------------------------------------------------------------
651 C-----------------------------------------------------------------------
652 C:::##INCLUDE '~/charmm_fcm/ctitla.fcm'
653 INTEGER MAXTIT
654 PARAMETER (MAXTIT=32)
655 INTEGER NTITLA,NTITLB
656 CHARACTER*80 TITLEA,TITLEB
657 COMMON /NTITLA/ NTITLA,NTITLB
658 COMMON /CTITLA/ TITLEA(MAXTIT),TITLEB(MAXTIT)
659 C..##IF SAVEFCM
660 C..##ENDIF
661 C-----------------------------------------------------------------------
662 C Passed variables
663 INTEGER NAT3,NADD,NPAR,NFREG,NFRET,BLATOM
664 INTEGER ATMPAR(2,*),ATMPAS(2,*),ATMPAD(2,*)
665 INTEGER BNBND(*),BIMAG(*)
666 INTEGER INBCMP(*),JNBCMP(*),PARDIM
667 INTEGER ITMX,IUNMOD,IUNRMD,SAVF
668 INTEGER NBOND,IB(*),JB(*)
669 REAL(KIND=8) X(*),Y(*),Z(*),AMASS(*),DDSCR(*)
670 REAL(KIND=8) DDV(NAT3,*),PARDDV(PARDIM,*),DDM(*),DDS(*)
671 REAL(KIND=8) DDF(*),PARDDF(*),DDEV(*),PARDDE(*)
672 REAL(KIND=8) DD1BLK(*),DD1BLL(*),DD1CMP(*)
673 REAL(KIND=8) TOLDIM,DDVALM
674 REAL(KIND=8) PARFRQ,CUTF1
675 LOGICAL LNOMA,LRAISE,LSCI,LBIG
676 C Local variables
677 INTEGER NATOM,NATP,NDIM,I,J,II,OLDFAS,OLDPRN,IUPD
678 INTEGER NPARC,NPARD,NPARS,NFCUT1,NFREG2,NFREG6
679 INTEGER IH1,IH2,IH3,IH4,IH5,IH6,IH7,IH8
680 INTEGER IS1,IS2,IS3,IS4,JSPACE,JSP,DDSS,DD5
681 INTEGER ISTRT,ISTOP,IPA1,IPA2,IRESF
682 INTEGER ATMPAF,INIDS,TRAROT
683 INTEGER SUBLIS,ATMCOR
684 INTEGER NFRRES,DDVBAS
685 INTEGER DDV2,DDVAL
686 INTEGER LENCM,NTR,NFRE,NFC,N1,N2,NFCUT,NSUBP
687 INTEGER SCIFV1,SCIFV2,SCIFV3,SCIFV4,SCIFV6
688 INTEGER DRATQ,ERATQ,E2RATQ,BDRATQ,INRATQ
689 INTEGER I620,I640,I660,I700,I720,I760,I800,I840,I880,I920
690 REAL(KIND=8) CVGMX,TOLER
691 LOGICAL LCARD,LAPPE,LPURG,LWDINI,QCALC,QMASWT,QMIX,QDIAG
692 C Begin
693 QCALC=.TRUE.
694 LWDINI=.FALSE.
695 INIDS=0
696 IS3=0
697 IS4=0
698 LPURG=.TRUE.
699 ITER=0
700 NADD=0
701 NFSAV=0
702 TOLER=TENM5
703 QDIAG=.TRUE.
704 CVGMX=HUNDRD
705 QMIX=.FALSE.
706 NATOM=NAT3/3
707 NFREG6=(NFREG-6)/NPAR
708 NFREG2=NFREG/2
709 NFRRES=(NFREG+6)/2
710 IF(NFREG.GT.PARDIM) CALL WRNDIE(-3,'<NMDIMB>',
711 1 'NFREG IS LARGER THAN PARDIM*3')
713 C ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
714 ASSIGN 801 TO I800 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
715 GOTO 800
716 801 CONTINUE
717 C ALLOCATE-SPACE-FOR-DIAGONALIZATION
718 ASSIGN 721 TO I720 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
719 GOTO 720
720 721 CONTINUE
721 C ALLOCATE-SPACE-FOR-REDUCED-BASIS
722 ASSIGN 761 TO I760 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
723 GOTO 760
724 761 CONTINUE
725 C ALLOCATE-SPACE-FOR-OTHER-ARRAYS
726 ASSIGN 921 TO I920 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
727 GOTO 920
728 921 CONTINUE
730 C Space allocation for working arrays of EISPACK
731 C diagonalization subroutines
732 IF(LSCI) THEN
733 C ALLOCATE-SPACE-FOR-LSCI
734 ASSIGN 841 TO I840 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
735 GOTO 840
736 841 CONTINUE
737 ELSE
738 C ALLOCATE-DUMMY-SPACE-FOR-LSCI
739 ASSIGN 881 TO I880 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
740 GOTO 880
741 881 CONTINUE
742 ENDIF
743 QMASWT=(.NOT.LNOMA)
744 IF(.NOT. QDISK) THEN
745 LENCM=INBCMP(NATOM-1)*9+NATOM*6
746 DO I=1,LENCM
747 DD1CMP(I)=0.0
748 ENDDO
749 OLDFAS=LFAST
750 QCMPCT=.TRUE.
751 LFAST = -1
752 CALL ENERGY(X,Y,Z,DX,DY,DZ,BNBND,BIMAG,NAT3,DD1CMP,.TRUE.,1)
753 LFAST=OLDFAS
754 QCMPCT=.FALSE.
756 C Mass weight DD1CMP matrix
758 CALL MASSDD(DD1CMP,DDM,INBCMP,JNBCMP,NATOM)
759 ELSE
760 CALL WRNDIE(-3,'<NMDIMB>','QDISK OPTION NOT SUPPORTED YET')
761 C DO I=1,LENDSK
762 C DD1CMP(I)=0.0
763 C ENDDO
764 C OLDFAS=LFAST
765 C LFAST = -1
766 ENDIF
768 C Fill DDV with six translation-rotation vectors
770 CALL TRROT(X,Y,Z,DDV,NAT3,1,DDM)
771 CALL CPARAY(HEAP(TRAROT),DDV,NAT3,1,6,1)
772 NTR=6
773 OLDPRN=PRNLEV
774 PRNLEV=1
775 CALL ORTHNM(1,6,NTR,HEAP(TRAROT),NAT3,.FALSE.,TOLER)
776 PRNLEV=OLDPRN
777 IF(IUNRMD .LT. 0) THEN
779 C If no previous basis is read
781 IF(PRNLEV.GE.2) WRITE(OUTU,502) NPAR
782 502 FORMAT(/' NMDIMB: Calculating initial basis from block ',
783 1 'diagonals'/' NMDIMB: The number of blocks is ',I5/)
784 NFRET = 6
785 DO I=1,NPAR
786 IS1=ATMPAR(1,I)
787 IS2=ATMPAR(2,I)
788 NDIM=(IS2-IS1+1)*3
789 NFRE=NDIM
790 IF(NFRE.GT.NFREG6) NFRE=NFREG6
791 IF(NFREG6.EQ.0) NFRE=1
792 CALL FILUPT(HEAP(IUPD),NDIM)
793 CALL MAKDDU(DD1BLK,DD1CMP,INBCMP,JNBCMP,HEAP(IUPD),
794 1 IS1,IS2,NATOM)
795 IF(PRNLEV.GE.9) CALL PRINTE(OUTU,EPROP,ETERM,'VIBR',
796 1 'ENR',.TRUE.,1,ZERO,ZERO)
798 C Generate the lower section of the matrix and diagonalize
800 C..##IF EISPACK
801 C..##ENDIF
802 IH1=1
803 NATP=NDIM+1
804 IH2=IH1+NATP
805 IH3=IH2+NATP
806 IH4=IH3+NATP
807 IH5=IH4+NATP
808 IH6=IH5+NATP
809 IH7=IH6+NATP
810 IH8=IH7+NATP
811 CALL DIAGQ(NDIM,NFRE,DD1BLK,PARDDV,DDS(IH2),DDS(IH3),
812 1 DDS(IH4),DDS(IH5),DDS,DDS(IH6),DDS(IH7),DDS(IH8),NADD)
813 C..##IF EISPACK
814 C..##ENDIF
816 C Put the PARDDV vectors into DDV and replace the elements which do
817 C not belong to the considered partitioned region by zeros.
819 CALL ADJNME(DDV,PARDDV,NAT3,NDIM,NFRE,NFRET,IS1,IS2)
820 IF(LSCI) THEN
821 DO J=1,NFRE
822 PARDDF(J)=CNVFRQ*SQRT(ABS(PARDDE(J)))
823 IF(PARDDE(J) .LT. 0.0) PARDDF(J)=-PARDDF(J)
824 ENDDO
825 ELSE
826 DO J=1,NFRE
827 PARDDE(J)=DDS(J)
828 PARDDF(J)=CNVFRQ*SQRT(ABS(PARDDE(J)))
829 IF(PARDDE(J) .LT. 0.0) PARDDF(J)=-PARDDF(J)
830 ENDDO
831 ENDIF
832 IF(PRNLEV.GE.2) THEN
833 WRITE(OUTU,512) I
834 WRITE(OUTU,514)
835 WRITE(OUTU,516) (J,PARDDF(J),J=1,NFRE)
836 ENDIF
837 NFRET=NFRET+NFRE
838 IF(NFRET .GE. NFREG) GOTO 10
839 ENDDO
840 512 FORMAT(/' NMDIMB: Diagonalization of part',I5,' completed')
841 514 FORMAT(' NMDIMB: Frequencies'/)
842 516 FORMAT(5(I4,F12.6))
843 10 CONTINUE
845 C Orthonormalize the eigenvectors
847 OLDPRN=PRNLEV
848 PRNLEV=1
849 CALL ORTHNM(1,NFRET,NFRET,DDV,NAT3,LPURG,TOLER)
850 PRNLEV=OLDPRN
852 C Do reduced basis diagonalization using the DDV vectors
853 C and get eigenvectors of zero iteration
855 IF(PRNLEV.GE.2) THEN
856 WRITE(OUTU,521) ITER
857 WRITE(OUTU,523) NFRET
858 ENDIF
859 521 FORMAT(/' NMDIMB: Iteration number = ',I5)
860 523 FORMAT(' NMDIMB: Dimension of the reduced basis set = ',I5)
861 IF(LBIG) THEN
862 IF(PRNLEV.GE.2) WRITE(OUTU,585) NFRET,IUNMOD
863 525 FORMAT(' NMDIMB: ',I5,' basis vectors are saved in unit',I5)
864 REWIND (UNIT=IUNMOD)
865 LCARD=.FALSE.
866 CALL WRTNMD(LCARD,1,NFRET,NAT3,DDV,DDSCR,DDEV,IUNMOD,AMASS)
867 CALL SAVEIT(IUNMOD)
868 ELSE
869 CALL CPARAY(HEAP(DDVBAS),DDV,NAT3,1,NFRET,1)
870 ENDIF
871 CALL RBDG(X,Y,Z,NAT3,NDIM,NFRET,DDV,DDF,DDEV,
872 1 DDSCR,HEAP(DD5),HEAP(DDSS),HEAP(DDV2),NADD,
873 2 INBCMP,JNBCMP,HEAP(DDVBAS),DD1CMP,QMIX,0,0,IS3,IS4,
874 3 CUTF1,NFCUT1,NFREG,HEAP(IUPD),DD1BLL,HEAP(SCIFV1),
875 4 HEAP(SCIFV2),HEAP(SCIFV3),HEAP(SCIFV4),HEAP(SCIFV6),
876 5 HEAP(DRATQ),HEAP(ERATQ),HEAP(E2RATQ),
877 6 HEAP(BDRATQ),HEAP(INRATQ),LSCI,LBIG,IUNMOD)
879 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
881 ASSIGN 621 TO I620 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
882 GOTO 620
883 621 CONTINUE
884 C SAVE-MODES
885 ASSIGN 701 TO I700 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
886 GOTO 700
887 701 CONTINUE
888 IF(ITER.EQ.ITMX) THEN
889 CALL CLEANHP(NAT3,NFREG,NPARD,NSUBP,PARDIM,DDV2,DDSS,DDVBAS,
890 1 DDVAL,JSPACE,TRAROT,
891 2 SCIFV1,SCIFV2,SCIFV3,SCIFV4,SCIFV6,
892 3 DRATQ,ERATQ,E2RATQ,BDRATQ,INRATQ,IUPD,ATMPAF,
893 4 ATMCOR,SUBLIS,LSCI,QDW,LBIG)
894 RETURN
895 ENDIF
896 ELSE
898 C Read in existing basis
900 IF(PRNLEV.GE.2) THEN
901 WRITE(OUTU,531)
902 531 FORMAT(/' NMDIMB: Calculations restarted')
903 ENDIF
904 C READ-MODES
905 ISTRT=1
906 ISTOP=99999999
907 LCARD=.FALSE.
908 LAPPE=.FALSE.
909 CALL RDNMD(LCARD,NFRET,NFREG,NAT3,NDIM,
910 1 DDV,DDSCR,DDF,DDEV,
911 2 IUNRMD,LAPPE,ISTRT,ISTOP)
912 NFRET=NDIM
913 IF(NFRET.GT.NFREG) THEN
914 NFRET=NFREG
915 CALL WRNDIE(-1,'<NMDIMB>',
916 1 'Not enough space to hold the basis. Increase NMODes')
917 ENDIF
918 C PRINT-MODES
919 IF(PRNLEV.GE.2) THEN
920 WRITE(OUTU,533) NFRET,IUNRMD
921 WRITE(OUTU,514)
922 WRITE(OUTU,516) (J,DDF(J),J=1,NFRET)
923 ENDIF
924 533 FORMAT(/' NMDIMB: ',I5,' restart modes read from unit ',I5)
925 NFRRES=NFRET
926 ENDIF
928 C -------------------------------------------------
929 C Here starts the mixed-basis diagonalization part.
930 C -------------------------------------------------
933 C Check cut-off frequency
935 CALL SELNMD(DDF,NFRET,CUTF1,NFCUT1)
936 C TEST-NFCUT1
937 IF(IUNRMD.LT.0) THEN
938 IF(NFCUT1*2-6.GT.NFREG) THEN
939 IF(PRNLEV.GE.2) WRITE(OUTU,537) DDF(NFRRES)
940 NFCUT1=NFRRES
941 CUTF1=DDF(NFRRES)
942 ENDIF
943 ELSE
944 CUTF1=DDF(NFRRES)
945 ENDIF
946 537 FORMAT(/' NMDIMB: Too many vectors for the given cutoff frequency'
947 1 /' Cutoff frequency is decreased to',F9.3)
949 C Compute the new partioning of the molecule
951 CALL PARTIC(NAT3,NFREG,NFCUT1,NPARMX,NPARC,ATMPAR,NFRRES,
952 1 PARDIM)
953 NPARS=NPARC
954 DO I=1,NPARC
955 ATMPAS(1,I)=ATMPAR(1,I)
956 ATMPAS(2,I)=ATMPAR(2,I)
957 ENDDO
958 IF(QDW) THEN
959 IF(IPAR1.EQ.0.OR.IPAR2.EQ.0) LWDINI=.TRUE.
960 IF(IPAR1.GE.IPAR2) LWDINI=.TRUE.
961 IF(IABS(IPAR1).GT.NPARC*2) LWDINI=.TRUE.
962 IF(IABS(IPAR2).GT.NPARC*2) LWDINI=.TRUE.
963 IF(ITER.EQ.0) LWDINI=.TRUE.
964 ENDIF
965 ITMX=ITMX+ITER
966 IF(PRNLEV.GE.2) THEN
967 WRITE(OUTU,543) ITER,ITMX
968 IF(QDW) WRITE(OUTU,545) IPAR1,IPAR2
969 ENDIF
970 543 FORMAT(/' NMDIMB: Previous iteration number = ',I8/
971 1 ' NMDIMB: Iteration number to reach = ',I8)
972 545 FORMAT(' NMDIMB: Previous sub-blocks = ',I5,2X,I5)
974 IF(SAVF.LE.0) SAVF=NPARC
975 IF(PRNLEV.GE.2) WRITE(OUTU,547) SAVF
976 547 FORMAT(' NMDIMB: Eigenvectors will be saved every',I5,
977 1 ' iterations')
979 C If double windowing is defined, the original block sizes are divided
980 C in two.
982 IF(QDW) THEN
983 NSUBP=1
984 CALL PARTID(NPARC,ATMPAR,NPARD,ATMPAD,NPARMX)
985 ATMPAF=ALLHP(INTEG4(NPARD*NPARD))
986 ATMCOR=ALLHP(INTEG4(NATOM))
987 DDVAL=ALLHP(IREAL8(NPARD*NPARD))
988 CALL CORARR(ATMPAD,NPARD,HEAP(ATMCOR),NATOM)
989 CALL PARLIS(HEAP(ATMCOR),HEAP(ATMPAF),INBCMP,JNBCMP,NPARD,
990 2 NSUBP,NATOM,X,Y,Z,NBOND,IB,JB,DD1CMP,HEAP(DDVAL),DDVALM)
991 SUBLIS=ALLHP(INTEG4(NSUBP*2))
992 CALL PARINT(HEAP(ATMPAF),NPARD,HEAP(SUBLIS),NSUBP)
993 CALL INIPAF(HEAP(ATMPAF),NPARD)
995 C Find out with which block to continue (double window method only)
997 IPA1=IPAR1
998 IPA2=IPAR2
999 IRESF=0
1000 IF(LWDINI) THEN
1001 ITER=0
1002 LWDINI=.FALSE.
1003 GOTO 500
1004 ENDIF
1005 DO II=1,NSUBP
1006 CALL IPART(HEAP(SUBLIS),II,IPAR1,IPAR2,HEAP(ATMPAF),
1007 1 NPARD,QCALC)
1008 IF((IPAR1.EQ.IPA1).AND.(IPAR2.EQ.IPA2)) GOTO 500
1009 ENDDO
1010 ENDIF
1011 500 CONTINUE
1013 C Main loop.
1015 DO WHILE((CVGMX.GT.TOLDIM).AND.(ITER.LT.ITMX))
1016 IF(.NOT.QDW) THEN
1017 ITER=ITER+1
1018 IF(PRNLEV.GE.2) WRITE(OUTU,553) ITER
1019 553 FORMAT(/' NMDIMB: Iteration number = ',I8)
1020 IF(INIDS.EQ.0) THEN
1021 INIDS=1
1022 ELSE
1023 INIDS=0
1024 ENDIF
1025 CALL PARTDS(NAT3,NPARC,ATMPAR,NPARS,ATMPAS,INIDS,NPARMX,
1026 1 DDF,NFREG,CUTF1,PARDIM,NFCUT1)
1027 C DO-THE-DIAGONALISATIONS
1028 ASSIGN 641 to I640 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1029 GOTO 640
1030 641 CONTINUE
1031 QDIAG=.FALSE.
1032 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1033 ASSIGN 622 TO I620 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1034 GOTO 620
1035 622 CONTINUE
1036 QDIAG=.TRUE.
1037 C SAVE-MODES
1038 ASSIGN 702 TO I700 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1039 GOTO 700
1040 702 CONTINUE
1042 ELSE
1043 DO II=1,NSUBP
1044 CALL IPART(HEAP(SUBLIS),II,IPAR1,IPAR2,HEAP(ATMPAF),
1045 1 NPARD,QCALC)
1046 IF(QCALC) THEN
1047 IRESF=IRESF+1
1048 ITER=ITER+1
1049 IF(PRNLEV.GE.2) WRITE(OUTU,553) ITER
1050 C DO-THE-DWIN-DIAGONALISATIONS
1051 ASSIGN 661 TO I660 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1052 GOTO 660
1053 661 CONTINUE
1054 ENDIF
1055 IF((IRESF.EQ.SAVF).OR.(ITER.EQ.ITMX)) THEN
1056 IRESF=0
1057 QDIAG=.FALSE.
1058 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1059 ASSIGN 623 TO I620 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1060 GOTO 620
1061 623 CONTINUE
1062 QDIAG=.TRUE.
1063 IF((CVGMX.LE.TOLDIM).OR.(ITER.EQ.ITMX)) GOTO 600
1064 C SAVE-MODES
1065 ASSIGN 703 TO I700 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1066 GOTO 700
1067 703 CONTINUE
1068 ENDIF
1069 ENDDO
1070 ENDIF
1071 ENDDO
1072 600 CONTINUE
1074 C SAVE-MODES
1075 ASSIGN 704 TO I700 ! { dg-warning "Obsolete: ASSIGN" "Obsolete: ASSIGN" }
1076 GOTO 700
1077 704 CONTINUE
1078 CALL CLEANHP(NAT3,NFREG,NPARD,NSUBP,PARDIM,DDV2,DDSS,DDVBAS,
1079 1 DDVAL,JSPACE,TRAROT,
1080 2 SCIFV1,SCIFV2,SCIFV3,SCIFV4,SCIFV6,
1081 3 DRATQ,ERATQ,E2RATQ,BDRATQ,INRATQ,IUPD,ATMPAF,
1082 4 ATMCOR,SUBLIS,LSCI,QDW,LBIG)
1083 RETURN
1084 C-----------------------------------------------------------------------
1085 C INTERNAL PROCEDURES
1086 C-----------------------------------------------------------------------
1087 C TO DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1088 620 CONTINUE
1089 IF(IUNRMD.LT.0) THEN
1090 CALL SELNMD(DDF,NFRET,CUTF1,NFC)
1091 N1=NFCUT1
1092 N2=(NFRET+6)/2
1093 NFCUT=MAX(N1,N2)
1094 IF(NFCUT*2-6 .GT. NFREG) THEN
1095 NFCUT=(NFREG+6)/2
1096 CUTF1=DDF(NFCUT)
1097 IF(PRNLEV.GE.2) THEN
1098 WRITE(OUTU,562) ITER
1099 WRITE(OUTU,564) CUTF1
1100 ENDIF
1101 ENDIF
1102 ELSE
1103 NFCUT=NFRET
1104 NFC=NFRET
1105 ENDIF
1106 562 FORMAT(/' NMDIMB: Not enough space to hold the residual vectors'/
1107 1 ' into DDV array during iteration ',I5)
1108 564 FORMAT(' Cutoff frequency is changed to ',F9.3)
1110 C do reduced diagonalization with preceding eigenvectors plus
1111 C residual vectors
1113 ISTRT=1
1114 ISTOP=NFCUT
1115 CALL CLETR(DDV,HEAP(TRAROT),NAT3,ISTRT,ISTOP,NFCUT,DDEV,DDF)
1116 CALL RNMTST(DDV,HEAP(DDVBAS),NAT3,DDSCR,DD1CMP,INBCMP,JNBCMP,
1117 2 7,NFCUT,CVGMX,NFCUT,NFC,QDIAG,LBIG,IUNMOD)
1118 NFSAV=NFCUT
1119 IF(QDIAG) THEN
1120 NFRET=NFCUT*2-6
1121 IF(PRNLEV.GE.2) WRITE(OUTU,566) NFRET
1122 566 FORMAT(/' NMDIMB: Diagonalization with residual vectors. '/
1123 1 ' Dimension of the reduced basis set'/
1124 2 ' before orthonormalization = ',I5)
1125 NFCUT=NFRET
1126 OLDPRN=PRNLEV
1127 PRNLEV=1
1128 CALL ORTHNM(1,NFRET,NFCUT,DDV,NAT3,LPURG,TOLER)
1129 PRNLEV=OLDPRN
1130 NFRET=NFCUT
1131 IF(PRNLEV.GE.2) WRITE(OUTU,568) NFRET
1132 568 FORMAT(' after orthonormalization = ',I5)
1133 IF(LBIG) THEN
1134 IF(PRNLEV.GE.2) WRITE(OUTU,570) NFCUT,IUNMOD
1135 570 FORMAT(' NMDIMB: ',I5,' basis vectors are saved in unit',I5)
1136 REWIND (UNIT=IUNMOD)
1137 LCARD=.FALSE.
1138 CALL WRTNMD(LCARD,1,NFCUT,NAT3,DDV,DDSCR,DDEV,IUNMOD,AMASS)
1139 CALL SAVEIT(IUNMOD)
1140 ELSE
1141 CALL CPARAY(HEAP(DDVBAS),DDV,NAT3,1,NFCUT,1)
1142 ENDIF
1143 QMIX=.FALSE.
1144 CALL RBDG(X,Y,Z,NAT3,NDIM,NFRET,DDV,DDF,DDEV,
1145 1 DDSCR,HEAP(DD5),HEAP(DDSS),HEAP(DDV2),NADD,
1146 2 INBCMP,JNBCMP,HEAP(DDVBAS),DD1CMP,QMIX,0,0,IS3,IS4,
1147 3 CUTF1,NFCUT1,NFREG,HEAP(IUPD),DD1BLL,HEAP(SCIFV1),
1148 4 HEAP(SCIFV2),HEAP(SCIFV3),HEAP(SCIFV4),HEAP(SCIFV6),
1149 5 HEAP(DRATQ),HEAP(ERATQ),HEAP(E2RATQ),
1150 6 HEAP(BDRATQ),HEAP(INRATQ),LSCI,LBIG,IUNMOD)
1151 CALL SELNMD(DDF,NFRET,CUTF1,NFCUT1)
1152 ENDIF
1153 GOTO I620 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1155 C-----------------------------------------------------------------------
1156 C TO DO-THE-DIAGONALISATIONS
1157 640 CONTINUE
1158 DO I=1,NPARC
1159 NFCUT1=NFRRES
1160 IS1=ATMPAR(1,I)
1161 IS2=ATMPAR(2,I)
1162 NDIM=(IS2-IS1+1)*3
1163 IF(PRNLEV.GE.2) WRITE(OUTU,573) I,IS1,IS2
1164 573 FORMAT(/' NMDIMB: Mixed diagonalization, part ',I5/
1165 1 ' NMDIMB: Block limits: ',I5,2X,I5)
1166 IF(NDIM+NFCUT1.GT.PARDIM) CALL WRNDIE(-3,'<NMDIMB>',
1167 1 'Error in dimension of block')
1168 NFRET=NFCUT1
1169 IF(NFRET.GT.NFREG) NFRET=NFREG
1170 CALL CLETR(DDV,HEAP(TRAROT),NAT3,1,NFCUT1,NFCUT,DDEV,DDF)
1171 NFCUT1=NFCUT
1172 CALL ADZER(DDV,1,NFCUT1,NAT3,IS1,IS2)
1173 NFSAV=NFCUT1
1174 OLDPRN=PRNLEV
1175 PRNLEV=1
1176 CALL ORTHNM(1,NFCUT1,NFCUT,DDV,NAT3,LPURG,TOLER)
1177 PRNLEV=OLDPRN
1178 CALL CPARAY(HEAP(DDVBAS),DDV,NAT3,1,NFCUT,1)
1179 NFRET=NDIM+NFCUT
1180 QMIX=.TRUE.
1181 CALL RBDG(X,Y,Z,NAT3,NDIM,NFRET,DDV,DDF,DDEV,
1182 1 DDSCR,HEAP(DD5),HEAP(DDSS),HEAP(DDV2),NADD,
1183 2 INBCMP,JNBCMP,HEAP(DDVBAS),DD1CMP,QMIX,IS1,IS2,IS3,IS4,
1184 3 CUTF1,NFCUT,NFREG,HEAP(IUPD),DD1BLL,HEAP(SCIFV1),
1185 4 HEAP(SCIFV2),HEAP(SCIFV3),HEAP(SCIFV4),HEAP(SCIFV6),
1186 5 HEAP(DRATQ),HEAP(ERATQ),HEAP(E2RATQ),
1187 6 HEAP(BDRATQ),HEAP(INRATQ),LSCI,LBIG,IUNMOD)
1188 QMIX=.FALSE.
1189 IF(NFCUT.GT.NFRRES) NFCUT=NFRRES
1190 NFCUT1=NFCUT
1191 NFRET=NFCUT
1192 ENDDO
1193 GOTO I640 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1195 C-----------------------------------------------------------------------
1196 C TO DO-THE-DWIN-DIAGONALISATIONS
1197 660 CONTINUE
1199 C Store the DDV vectors into DDVBAS
1201 NFCUT1=NFRRES
1202 IS1=ATMPAD(1,IPAR1)
1203 IS2=ATMPAD(2,IPAR1)
1204 IS3=ATMPAD(1,IPAR2)
1205 IS4=ATMPAD(2,IPAR2)
1206 NDIM=(IS2-IS1+IS4-IS3+2)*3
1207 IF(PRNLEV.GE.2) WRITE(OUTU,577) IPAR1,IPAR2,IS1,IS2,IS3,IS4
1208 577 FORMAT(/' NMDIMB: Mixed double window diagonalization, parts ',
1209 1 2I5/
1210 2 ' NMDIMB: Block limits: ',I5,2X,I5,4X,I5,2X,I5)
1211 IF(NDIM+NFCUT1.GT.PARDIM) CALL WRNDIE(-3,'<NMDIMB>',
1212 1 'Error in dimension of block')
1213 NFRET=NFCUT1
1214 IF(NFRET.GT.NFREG) NFRET=NFREG
1216 C Prepare the DDV vectors consisting of 6 translations-rotations
1217 C + eigenvectors from 7 to NFCUT1 + cartesian displacements vectors
1218 C spanning the atoms from IS1 to IS2
1220 CALL CLETR(DDV,HEAP(TRAROT),NAT3,1,NFCUT1,NFCUT,DDEV,DDF)
1221 NFCUT1=NFCUT
1222 NFSAV=NFCUT1
1223 CALL ADZERD(DDV,1,NFCUT1,NAT3,IS1,IS2,IS3,IS4)
1224 OLDPRN=PRNLEV
1225 PRNLEV=1
1226 CALL ORTHNM(1,NFCUT1,NFCUT,DDV,NAT3,LPURG,TOLER)
1227 PRNLEV=OLDPRN
1228 CALL CPARAY(HEAP(DDVBAS),DDV,NAT3,1,NFCUT,1)
1230 NFRET=NDIM+NFCUT
1231 QMIX=.TRUE.
1232 CALL RBDG(X,Y,Z,NAT3,NDIM,NFRET,DDV,DDF,DDEV,
1233 1 DDSCR,HEAP(DD5),HEAP(DDSS),HEAP(DDV2),NADD,
1234 2 INBCMP,JNBCMP,HEAP(DDVBAS),DD1CMP,QMIX,IS1,IS2,IS3,IS4,
1235 3 CUTF1,NFCUT,NFREG,HEAP(IUPD),DD1BLL,HEAP(SCIFV1),
1236 4 HEAP(SCIFV2),HEAP(SCIFV3),HEAP(SCIFV4),HEAP(SCIFV6),
1237 5 HEAP(DRATQ),HEAP(ERATQ),HEAP(E2RATQ),
1238 6 HEAP(BDRATQ),HEAP(INRATQ),LSCI,LBIG,IUNMOD)
1239 QMIX=.FALSE.
1241 IF(NFCUT.GT.NFRRES) NFCUT=NFRRES
1242 NFCUT1=NFCUT
1243 NFRET=NFCUT
1244 GOTO I660 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1246 C-----------------------------------------------------------------------
1247 C TO SAVE-MODES
1248 700 CONTINUE
1249 IF(PRNLEV.GE.2) WRITE(OUTU,583) IUNMOD
1250 583 FORMAT(/' NMDIMB: Saving the eigenvalues and eigenvectors to unit'
1251 1 ,I4)
1252 REWIND (UNIT=IUNMOD)
1253 ISTRT=1
1254 ISTOP=NFSAV
1255 LCARD=.FALSE.
1256 IF(PRNLEV.GE.2) WRITE(OUTU,585) NFSAV,IUNMOD
1257 585 FORMAT(' NMDIMB: ',I5,' modes are saved in unit',I5)
1258 CALL WRTNMD(LCARD,ISTRT,ISTOP,NAT3,DDV,DDSCR,DDEV,IUNMOD,
1259 1 AMASS)
1260 CALL SAVEIT(IUNMOD)
1261 GOTO I700 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1263 C-----------------------------------------------------------------------
1264 C TO ALLOCATE-SPACE-FOR-DIAGONALIZATION
1265 720 CONTINUE
1266 DDV2=ALLHP(IREAL8((PARDIM+3)*(PARDIM+3)))
1267 JSPACE=IREAL8((PARDIM+4))*8
1268 JSP=IREAL8(((PARDIM+3)*(PARDIM+4))/2)
1269 JSPACE=JSPACE+JSP
1270 DDSS=ALLHP(JSPACE)
1271 DD5=DDSS+JSPACE-JSP
1272 GOTO I720 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1274 C-----------------------------------------------------------------------
1275 C TO ALLOCATE-SPACE-FOR-REDUCED-BASIS
1276 760 CONTINUE
1277 IF(LBIG) THEN
1278 DDVBAS=ALLHP(IREAL8(NAT3))
1279 ELSE
1280 DDVBAS=ALLHP(IREAL8(NFREG*NAT3))
1281 ENDIF
1282 GOTO I760 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1284 C-----------------------------------------------------------------------
1285 C TO ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
1286 800 CONTINUE
1287 TRAROT=ALLHP(IREAL8(6*NAT3))
1288 GOTO I800 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1290 C-----------------------------------------------------------------------
1291 C TO ALLOCATE-SPACE-FOR-LSCI
1292 840 CONTINUE
1293 SCIFV1=ALLHP(IREAL8(PARDIM+3))
1294 SCIFV2=ALLHP(IREAL8(PARDIM+3))
1295 SCIFV3=ALLHP(IREAL8(PARDIM+3))
1296 SCIFV4=ALLHP(IREAL8(PARDIM+3))
1297 SCIFV6=ALLHP(IREAL8(PARDIM+3))
1298 DRATQ=ALLHP(IREAL8(PARDIM+3))
1299 ERATQ=ALLHP(IREAL8(PARDIM+3))
1300 E2RATQ=ALLHP(IREAL8(PARDIM+3))
1301 BDRATQ=ALLHP(IREAL8(PARDIM+3))
1302 INRATQ=ALLHP(INTEG4(PARDIM+3))
1303 GOTO I840 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1305 C-----------------------------------------------------------------------
1306 C TO ALLOCATE-DUMMY-SPACE-FOR-LSCI
1307 880 CONTINUE
1308 SCIFV1=ALLHP(IREAL8(2))
1309 SCIFV2=ALLHP(IREAL8(2))
1310 SCIFV3=ALLHP(IREAL8(2))
1311 SCIFV4=ALLHP(IREAL8(2))
1312 SCIFV6=ALLHP(IREAL8(2))
1313 DRATQ=ALLHP(IREAL8(2))
1314 ERATQ=ALLHP(IREAL8(2))
1315 E2RATQ=ALLHP(IREAL8(2))
1316 BDRATQ=ALLHP(IREAL8(2))
1317 INRATQ=ALLHP(INTEG4(2))
1318 GOTO I880 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1320 C-----------------------------------------------------------------------
1321 C TO ALLOCATE-SPACE-FOR-OTHER-ARRAYS
1322 920 CONTINUE
1323 IUPD=ALLHP(INTEG4(PARDIM+3))
1324 GOTO I920 ! { dg-warning "Obsolete: Assigned" "Assigned GO TO" }
1325 C.##ELSE
1326 C.##ENDIF