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