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