1c { dg-do compile }
2CHARMM Element source/dimb/nmdimb.src 1.1
3C.##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)
10C-----------------------------------------------------------------------
11C     01-Jul-1992 David Perahia, Liliane Mouawad
12C     15-Dec-1994 Herman van Vlijmen
13C
14C     This is the main routine for the mixed-basis diagonalization.
15C     See: L.Mouawad and D.Perahia, Biopolymers (1993), 33, 599,
16C     and: D.Perahia and L.Mouawad, Comput. Chem. (1995), 19, 241.
17C     The method iteratively solves the diagonalization of the
18C     Hessian matrix. To save memory space, it uses a compressed
19C     form of the Hessian, which only contains the nonzero elements.
20C     In the diagonalization process, approximate eigenvectors are
21C     mixed with Cartesian coordinates to form a reduced basis. The
22C     Hessian is then diagonalized in the reduced basis. By iterating
23C     over different sets of Cartesian coordinates the method ultimately
24C     converges to the exact eigenvalues and eigenvectors (up to the
25C     requested accuracy).
26C     If no existing basis set is read, an initial basis will be created
27C     which consists of the low-frequency eigenvectors of diagonal blocks
28C     of the Hessian.
29C-----------------------------------------------------------------------
30C-----------------------------------------------------------------------
31C:::##INCLUDE '~/charmm_fcm/impnon.fcm'
32C..##IF VAX IRIS HPUX IRIS GNU CSPP OS2 GWS CRAY ALPHA
33      IMPLICIT NONE
34C..##ENDIF
35C-----------------------------------------------------------------------
36C-----------------------------------------------------------------------
37C:::##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
44C..##IF SAVEFCM
45C..##ENDIF
46C-----------------------------------------------------------------------
47C-----------------------------------------------------------------------
48C:::##INCLUDE '~/charmm_fcm/dimens.fcm'
49      INTEGER LARGE,MEDIUM,SMALL,REDUCE
50C..##IF QUANTA
51C..##ELIF T3D
52C..##ELSE
53      PARAMETER (LARGE=60120, MEDIUM=25140, SMALL=6120)
54C..##ENDIF
55      PARAMETER (REDUCE=15000)
56      INTEGER SIZE
57C..##IF XLARGE
58C..##ELIF XXLARGE
59C..##ELIF LARGE
60C..##ELIF MEDIUM
61      PARAMETER (SIZE=MEDIUM)
62C..##ELIF REDUCE
63C..##ELIF SMALL
64C..##ELIF XSMALL
65C..##ENDIF
66C..##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)
75C..##ENDIF
76      INTEGER MAXCSP, MAXHSET
77C..##IF HMCM
78      PARAMETER (MAXHSET = 200)
79C..##ELSE
80C..##ENDIF
81C..##IF REDUCE
82C..##ELSE
83      PARAMETER (MAXCSP = 500)
84C..##ENDIF
85C..##IF HMCM
86      INTEGER MAXHCM,MAXPCM,MAXRCM
87C...##IF REDUCE
88C...##ELSE
89      PARAMETER (MAXHCM=500)
90      PARAMETER (MAXPCM=5000)
91      PARAMETER (MAXRCM=2000)
92C...##ENDIF
93C..##ENDIF
94      INTEGER MXCMSZ
95C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
96C..##ELSE
97      PARAMETER (MXCMSZ = 5000)
98C..##ENDIF
99      INTEGER CHRSIZ
100      PARAMETER (CHRSIZ = SIZE)
101      INTEGER MAXATB
102C..##IF REDUCE
103C..##ELIF QUANTA
104C..##ELSE
105      PARAMETER (MAXATB = 200)
106C..##ENDIF
107      INTEGER MAXVEC
108C..##IFN VECTOR PARVECT
109      PARAMETER (MAXVEC = 10)
110C..##ELIF LARGE XLARGE XXLARGE
111C..##ELIF MEDIUM
112C..##ELIF SMALL REDUCE
113C..##ELIF XSMALL
114C..##ELSE
115C..##ENDIF
116      INTEGER IATBMX
117      PARAMETER (IATBMX = 8)
118      INTEGER MAXHB
119C..##IF LARGE XLARGE XXLARGE
120C..##ELIF MEDIUM
121      PARAMETER (MAXHB = 8000)
122C..##ELIF SMALL
123C..##ELIF REDUCE XSMALL
124C..##ELSE
125C..##ENDIF
126      INTEGER MAXTRN,MAXSYM
127C..##IFN NOIMAGES
128      PARAMETER (MAXTRN = 5000)
129      PARAMETER (MAXSYM = 192)
130C..##ELSE
131C..##ENDIF
132C..##IF LONEPAIR (lonepair_max)
133      INTEGER MAXLP,MAXLPH
134C...##IF REDUCE
135C...##ELSE
136      PARAMETER (MAXLP  = 2000)
137      PARAMETER (MAXLPH = 4000)
138C...##ENDIF
139C..##ENDIF (lonepair_max)
140      INTEGER NOEMAX,NOEMX2
141C..##IF REDUCE
142C..##ELSE
143      PARAMETER (NOEMAX = 2000)
144      PARAMETER (NOEMX2 = 4000)
145C..##ENDIF
146      INTEGER MAXATC, MAXCB, MAXCH, MAXCI, MAXCP, MAXCT, MAXITC, MAXNBF
147C..##IF REDUCE
148C..##ELIF MMFF CFF
149      PARAMETER (MAXATC = 500, MAXCB = 1500, MAXCH = 3200, MAXCI = 600,
150     &           MAXCP  = 3000,MAXCT = 15500,MAXITC = 200, MAXNBF=1000)
151C..##ELIF YAMMP
152C..##ELIF LARGE
153C..##ELSE
154C..##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
160C..##IF LARGE XLARGE XXLARGE
161C..##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)
166C...##IF MCSS
167C...##ELSE
168      PARAMETER (MAXSEG = 1000)
169C...##ENDIF
170C..##ELIF SMALL
171C..##ELIF XSMALL
172C..##ELIF REDUCE
173C..##ELSE
174C..##ENDIF
175C..##IF NOIMAGES
176C..##ELSE
177      PARAMETER (MAXAIM = 2*SIZE)
178      PARAMETER (MAXGRP = 2*SIZE/3)
179C..##ENDIF
180      INTEGER REDMAX,REDMX2
181C..##IF REDUCE
182C..##ELSE
183      PARAMETER (REDMAX = 20)
184      PARAMETER (REDMX2 = 80)
185C..##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,
190C..##IF YAMMP
191C..##ELSE
192     &           MXRTX = 5000, MXRTHA = 300, MXRTHD = 300,
193C..##ENDIF
194     &           MXRTBL = 5000, NICM = 10)
195      INTEGER NMFTAB,  NMCTAB,  NMCATM,  NSPLIN
196C..##IF REDUCE
197C..##ELSE
198      PARAMETER (NMFTAB = 200, NMCTAB = 3, NMCATM = 12000, NSPLIN = 3)
199C..##ENDIF
200      INTEGER MAXSHK
201C..##IF XSMALL
202C..##ELIF REDUCE
203C..##ELSE
204      PARAMETER (MAXSHK = SIZE*3/4)
205C..##ENDIF
206      INTEGER SCRMAX
207C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
208C..##ELSE
209      PARAMETER (SCRMAX = 5000)
210C..##ENDIF
211C..##IF TSM
212      INTEGER MXPIGG
213C...##IF REDUCE
214C...##ELSE
215      PARAMETER (MXPIGG=500)
216C...##ENDIF
217      INTEGER MXCOLO,MXPUMB
218      PARAMETER (MXCOLO=20,MXPUMB=20)
219C..##ENDIF
220C..##IF ADUMB
221      INTEGER MAXUMP, MAXEPA, MAXNUM
222C...##IF REDUCE
223C...##ELSE
224      PARAMETER (MAXUMP = 10, MAXNUM = 4)
225C...##ENDIF
226C..##ENDIF
227      INTEGER MAXING
228      PARAMETER (MAXING=1000)
229C..##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)
236C..##ENDIF
237C-----------------------------------------------------------------------
238C-----------------------------------------------------------------------
239C:::##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
243C..##IF SINGLE
244C..##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)
251C..##ENDIF
252      REAL(KIND=8)     FIFTY, SIXTY, SVNTY2, EIGHTY, NINETY, HUNDRD,
253     &           ONE2TY, ONE8TY, THRHUN, THR6TY, NINE99, FIFHUN, THOSND,
254     &           FTHSND,MEGA
255C..##IF SINGLE
256C..##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)
262C..##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
268C..##IF SINGLE
269C..##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)
278C..##ENDIF
279      REAL(KIND=8) ANUM,FMARK
280      REAL(KIND=8) RSMALL,RBIG
281C..##IF SINGLE
282C..##ELSE
283      PARAMETER (ANUM=9999.0D0, FMARK=-999.0D0)
284      PARAMETER (RSMALL=1.0D-10,RBIG=1.0D20)
285C..##ENDIF
286      REAL(KIND=8) RPRECI,RBIGST
287C..##IF VAX DEC
288C..##ELIF IBM
289C..##ELIF CRAY
290C..##ELIF ALPHA T3D T3E
291C..##ELSE
292C...##IF SINGLE
293C...##ELSE
294      PARAMETER (RPRECI = 2.22045D-16, RBIGST = 4.49423D+307)
295C...##ENDIF
296C..##ENDIF
297C-----------------------------------------------------------------------
298C-----------------------------------------------------------------------
299C:::##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
311C..##IF AMBER
312C..##ELIF DISCOVER
313C..##ELSE
314      PARAMETER (CCELEC=332.0716D0)
315C..##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 )
328C..##IF MMFF
329      REAL(KIND=8) MDAKCAL
330      parameter(MDAKCAL=143.9325D0)
331C..##ENDIF
332      REAL(KIND=8) DEBYEC
333      PARAMETER ( DEBYEC = 2.541766D0 / BOHRR )
334      REAL(KIND=8) ZEROC
335      PARAMETER ( ZEROC = 298.15D0 )
336C-----------------------------------------------------------------------
337C-----------------------------------------------------------------------
338C:::##INCLUDE '~/charmm_fcm/exfunc.fcm'
339C..##IF ACE
340C..##ENDIF
341C..##IF ADUMB
342C..##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
354C..##IF ACE
355     *           ,GETNNB
356C..##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
362C..##IF ADUMB
363     *           ,UMFI
364C..##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
378C..##IF ADUMB
379     *           ,UMFI
380C..##ENDIF
381C..##IF ACE
382     *           ,GETNNB
383C..##ENDIF
384C..##IFN NOIMAGES
385      INTEGER IMATOM
386      EXTERNAL IMATOM
387C..##ENDIF
388C..##IF MBOND
389C..##ENDIF
390C..##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
408C..##ENDIF
409C-----------------------------------------------------------------------
410C-----------------------------------------------------------------------
411C:::##INCLUDE '~/charmm_fcm/stack.fcm'
412      INTEGER STKSIZ
413C..##IFN UNICOS
414C...##IF LARGE XLARGE
415C...##ELIF MEDIUM REDUCE
416      PARAMETER (STKSIZ=4000000)
417C...##ELIF SMALL
418C...##ELIF XSMALL
419C...##ELIF XXLARGE
420C...##ELSE
421C...##ENDIF
422      INTEGER LSTUSD,MAXUSD,STACK
423      COMMON /ISTACK/ LSTUSD,MAXUSD,STACK(STKSIZ)
424C..##ELSE
425C..##ENDIF
426C..##IF SAVEFCM
427C..##ENDIF
428C-----------------------------------------------------------------------
429C-----------------------------------------------------------------------
430C:::##INCLUDE '~/charmm_fcm/heap.fcm'
431      INTEGER HEAPDM
432C..##IFN UNICOS (unicos)
433C...##IF XXLARGE (size)
434C...##ELIF LARGE XLARGE (size)
435C...##ELIF MEDIUM (size)
436C....##IF T3D (t3d2)
437C....##ELIF TERRA (t3d2)
438C....##ELIF ALPHA (t3d2)
439C....##ELIF T3E (t3d2)
440C....##ELSE (t3d2)
441      PARAMETER (HEAPDM=2048000)
442C....##ENDIF (t3d2)
443C...##ELIF SMALL (size)
444C...##ELIF REDUCE (size)
445C...##ELIF XSMALL (size)
446C...##ELSE (size)
447C...##ENDIF (size)
448      INTEGER FREEHP,HEAPSZ,HEAP
449      COMMON /HEAPST/ FREEHP,HEAPSZ,HEAP(HEAPDM)
450      LOGICAL LHEAP(HEAPDM)
451      EQUIVALENCE (LHEAP,HEAP)
452C..##ELSE (unicos)
453C..##ENDIF (unicos)
454C..##IF SAVEFCM (save)
455C..##ENDIF (save)
456C-----------------------------------------------------------------------
457C-----------------------------------------------------------------------
458C:::##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
466C..##IF SAVEFCM
467C..##ENDIF
468C-----------------------------------------------------------------------
469C-----------------------------------------------------------------------
470C:::##INCLUDE '~/charmm_fcm/deriv.fcm'
471      REAL(KIND=8) DX,DY,DZ
472      COMMON /DERIVR/ DX(MAXAIM),DY(MAXAIM),DZ(MAXAIM)
473C..##IF SAVEFCM
474C..##ENDIF
475C-----------------------------------------------------------------------
476C-----------------------------------------------------------------------
477C:::##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
487C..##IF ACE
488     &      , SELF, SCREEN, COUL ,SOLV, INTER
489C..##ENDIF
490C..##IF FLUCQ
491     &       ,FQKIN
492C..##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
504C..##IF ACE
505     &         , SELF   = 45, SCREEN = 46, COUL   = 47,
506     &           SOLV   = 48, INTER  = 49
507C..##ENDIF
508C..##IF FLUCQ
509     &          ,FQKIN  = 50
510C..##ENDIF
511     &          )
512C..##IF ACE
513C..##ENDIF
514C..##IF GRID
515C..##ENDIF
516C..##IF FLUCQ
517C..##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
526C..##IF HMCM
527     &       , HMCM
528C..##ENDIF
529C..##IF ADUMB
530     &       , ADUMB
531C..##ENDIF
532     &       , HYDR
533C..##IF FLUCQ
534     &       , FQPOL
535C..##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
552C..##IF HMCM
553     &         , HMCM   = 61
554C..##ENDIF
555C..##IF ADUMB
556     &         , ADUMB  = 62
557C..##ENDIF
558     &         , HYDR   = 63
559C..##IF FLUCQ
560     &         , FQPOL  = 65
561C..##ENDIF
562     &           )
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)
585C..##IF SAVEFCM
586C..##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)
594C..##IF SAVEFCM
595C..##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
602C..##IF SAVEFCM
603C..##ENDIF
604C..##IF ACE
605C..##ENDIF
606C..##IF FLUCQ
607C..##ENDIF
608C..##IF ADUMB
609C..##ENDIF
610C..##IF GRID
611C..##ENDIF
612C..##IF FLUCQ
613C..##ENDIF
614C..##IF TSM
615      REAL(KIND=8) TSMTRM(LENENT),TSMTMP(LENENT)
616      COMMON /TSMENG/ TSMTRM,TSMTMP
617C...##IF SAVEFCM
618C...##ENDIF
619C..##ENDIF
620      REAL(KIND=8) EHQBM
621      LOGICAL HQBM
622      COMMON /HQBMVAR/HQBM
623C..##IF SAVEFCM
624C..##ENDIF
625C-----------------------------------------------------------------------
626C-----------------------------------------------------------------------
627C:::##INCLUDE '~/charmm_fcm/dimb.fcm'
628C..##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
647C...##IF SAVEFCM
648C...##ENDIF
649C..##ENDIF (dimbfcm)
650C-----------------------------------------------------------------------
651C-----------------------------------------------------------------------
652C:::##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)
659C..##IF SAVEFCM
660C..##ENDIF
661C-----------------------------------------------------------------------
662C 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
676C 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
692C 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')
712C
713C ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
714      ASSIGN 801 TO I800 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
715      GOTO 800
716 801  CONTINUE
717C ALLOCATE-SPACE-FOR-DIAGONALIZATION
718      ASSIGN 721 TO I720 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
719      GOTO 720
720 721  CONTINUE
721C ALLOCATE-SPACE-FOR-REDUCED-BASIS
722      ASSIGN 761 TO I760 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
723      GOTO 760
724 761  CONTINUE
725C ALLOCATE-SPACE-FOR-OTHER-ARRAYS
726      ASSIGN 921 TO I920 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
727      GOTO 920
728 921  CONTINUE
729C
730C Space allocation for working arrays of EISPACK
731C diagonalization subroutines
732      IF(LSCI) THEN
733C ALLOCATE-SPACE-FOR-LSCI
734         ASSIGN 841 TO I840 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
735         GOTO 840
736 841     CONTINUE
737      ELSE
738C ALLOCATE-DUMMY-SPACE-FOR-LSCI
739         ASSIGN 881 TO I880 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: 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.
755C
756C Mass weight DD1CMP matrix
757C
758         CALL MASSDD(DD1CMP,DDM,INBCMP,JNBCMP,NATOM)
759      ELSE
760         CALL WRNDIE(-3,'<NMDIMB>','QDISK OPTION NOT SUPPORTED YET')
761C         DO I=1,LENDSK
762C            DD1CMP(I)=0.0
763C         ENDDO
764C         OLDFAS=LFAST
765C         LFAST = -1
766      ENDIF
767C
768C Fill DDV with six translation-rotation vectors
769C
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
778C
779C If no previous basis is read
780C
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)
797C
798C Generate the lower section of the matrix and diagonalize
799C
800C..##IF EISPACK
801C..##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)
813C..##IF EISPACK
814C..##ENDIF
815C
816C Put the PARDDV vectors into DDV and replace the elements which do
817C not belong to the considered partitioned region by zeros.
818C
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
844C
845C Orthonormalize the eigenvectors
846C
847         OLDPRN=PRNLEV
848         PRNLEV=1
849         CALL ORTHNM(1,NFRET,NFRET,DDV,NAT3,LPURG,TOLER)
850         PRNLEV=OLDPRN
851C
852C Do reduced basis diagonalization using the DDV vectors
853C and get eigenvectors of zero iteration
854C
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)
878C
879C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
880C
881         ASSIGN 621 TO I620 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
882         GOTO 620
883 621     CONTINUE
884C SAVE-MODES
885         ASSIGN 701 TO I700 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: 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
897C
898C Read in existing basis
899C
900         IF(PRNLEV.GE.2) THEN
901            WRITE(OUTU,531)
902 531        FORMAT(/' NMDIMB: Calculations restarted')
903         ENDIF
904C 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
918C 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
927C
928C -------------------------------------------------
929C Here starts the mixed-basis diagonalization part.
930C -------------------------------------------------
931C
932C
933C Check cut-off frequency
934C
935      CALL SELNMD(DDF,NFRET,CUTF1,NFCUT1)
936C 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)
948C
949C Compute the new partioning of the molecule
950C
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)
973C
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')
978C
979C If double windowing is defined, the original block sizes are divided
980C in two.
981C
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)
994C
995C Find out with which block to continue (double window method only)
996C
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
1012C
1013C Main loop.
1014C
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)
1027C DO-THE-DIAGONALISATIONS
1028            ASSIGN 641 to I640 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
1029            GOTO 640
1030 641        CONTINUE
1031            QDIAG=.FALSE.
1032C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1033            ASSIGN 622 TO I620 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
1034            GOTO 620
1035 622        CONTINUE
1036            QDIAG=.TRUE.
1037C SAVE-MODES
1038            ASSIGN 702 TO I700 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
1039            GOTO 700
1040 702        CONTINUE
1041C
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
1050C DO-THE-DWIN-DIAGONALISATIONS
1051                  ASSIGN 661 TO I660 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: 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.
1058C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1059                  ASSIGN 623 TO I620 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
1060                  GOTO 620
1061 623              CONTINUE
1062                  QDIAG=.TRUE.
1063                  IF((CVGMX.LE.TOLDIM).OR.(ITER.EQ.ITMX)) GOTO 600
1064C SAVE-MODES
1065                  ASSIGN 703 TO I700 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: ASSIGN" }
1066                  GOTO 700
1067 703              CONTINUE
1068               ENDIF
1069            ENDDO
1070         ENDIF
1071      ENDDO
1072 600  CONTINUE
1073C
1074C SAVE-MODES
1075      ASSIGN 704 TO I700 ! { dg-warning "Deleted feature: ASSIGN" "Deleted feature: 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
1084C-----------------------------------------------------------------------
1085C INTERNAL PROCEDURES
1086C-----------------------------------------------------------------------
1087C 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)
1109C
1110C do reduced diagonalization with preceding eigenvectors plus
1111C residual vectors
1112C
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 "Deleted feature: Assigned" "Assigned GO TO" }
1154C
1155C-----------------------------------------------------------------------
1156C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1194C
1195C-----------------------------------------------------------------------
1196C TO DO-THE-DWIN-DIAGONALISATIONS
1197 660  CONTINUE
1198C
1199C Store the DDV vectors into DDVBAS
1200C
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
1215C
1216C Prepare the DDV vectors consisting of 6 translations-rotations
1217C + eigenvectors from 7 to NFCUT1 + cartesian displacements vectors
1218C spanning the atoms from IS1 to IS2
1219C
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)
1229C
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.
1240C
1241      IF(NFCUT.GT.NFRRES) NFCUT=NFRRES
1242      NFCUT1=NFCUT
1243      NFRET=NFCUT
1244      GOTO I660 ! { dg-warning "Deleted feature: Assigned" "Assigned GO TO" }
1245C
1246C-----------------------------------------------------------------------
1247C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1262C
1263C-----------------------------------------------------------------------
1264C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1273C
1274C-----------------------------------------------------------------------
1275C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1283C
1284C-----------------------------------------------------------------------
1285C TO ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
1286 800  CONTINUE
1287      TRAROT=ALLHP(IREAL8(6*NAT3))
1288      GOTO I800 ! { dg-warning "Deleted feature: Assigned" "Assigned GO TO" }
1289C
1290C-----------------------------------------------------------------------
1291C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1304C
1305C-----------------------------------------------------------------------
1306C 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 "Deleted feature: Assigned" "Assigned GO TO" }
1319C
1320C-----------------------------------------------------------------------
1321C TO ALLOCATE-SPACE-FOR-OTHER-ARRAYS
1322 920  CONTINUE
1323      IUPD=ALLHP(INTEG4(PARDIM+3))
1324      GOTO I920 ! { dg-warning "Deleted feature: Assigned" "Assigned GO TO" }
1325C.##ELSE
1326C.##ENDIF
1327      END
1328