1 CHARMM Element source/dimb/nmdimb.src 1.1
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
)
9 C-----------------------------------------------------------------------
10 C 01-Jul-1992 David Perahia, Liliane Mouawad
11 C 15-Dec-1994 Herman van Vlijmen
13 C This is the main routine for the mixed-basis diagonalization.
14 C See: L.Mouawad and D.Perahia, Biopolymers (1993), 33, 599,
15 C and: D.Perahia and L.Mouawad, Comput. Chem. (1995), 19, 241.
16 C The method iteratively solves the diagonalization of the
17 C Hessian matrix. To save memory space, it uses a compressed
18 C form of the Hessian, which only contains the nonzero elements.
19 C In the diagonalization process, approximate eigenvectors are
20 C mixed with Cartesian coordinates to form a reduced basis. The
21 C Hessian is then diagonalized in the reduced basis. By iterating
22 C over different sets of Cartesian coordinates the method ultimately
23 C converges to the exact eigenvalues and eigenvectors (up to the
24 C requested accuracy).
25 C If no existing basis set is read, an initial basis will be created
26 C which consists of the low-frequency eigenvectors of diagonal blocks
28 C-----------------------------------------------------------------------
29 C-----------------------------------------------------------------------
30 C:::##INCLUDE '~/charmm_fcm/impnon.fcm'
31 C..##IF VAX IRIS HPUX IRIS GNU CSPP OS2 GWS CRAY ALPHA
34 C-----------------------------------------------------------------------
35 C-----------------------------------------------------------------------
36 C:::##INCLUDE '~/charmm_fcm/stream.fcm'
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
45 C-----------------------------------------------------------------------
46 C-----------------------------------------------------------------------
47 C:::##INCLUDE '~/charmm_fcm/dimens.fcm'
48 INTEGER LARGE
,MEDIUM
,SMALL
,REDUCE
52 PARAMETER (LARGE
=60120, MEDIUM
=25140, SMALL
=6120)
54 PARAMETER (REDUCE
=15000)
60 PARAMETER (SIZE
=MEDIUM
)
67 parameter(MAXDEFI
=250)
68 INTEGER NAME0
,NAMEQ0
,NRES0
,KRES0
69 PARAMETER (NAME0
=4,NAMEQ0
=10,NRES0
=4,KRES0
=4)
73 PARAMETER (MAXAUX
= 10)
75 INTEGER MAXCSP
, MAXHSET
77 PARAMETER (MAXHSET
= 200)
82 PARAMETER (MAXCSP
= 500)
85 INTEGER MAXHCM
,MAXPCM
,MAXRCM
88 PARAMETER (MAXHCM
=500)
89 PARAMETER (MAXPCM
=5000)
90 PARAMETER (MAXRCM
=2000)
94 C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
96 PARAMETER (MXCMSZ
= 5000)
99 PARAMETER (CHRSIZ
= SIZE
)
104 PARAMETER (MAXATB
= 200)
107 C..##IFN VECTOR PARVECT
108 PARAMETER (MAXVEC
= 10)
109 C..##ELIF LARGE XLARGE XXLARGE
111 C..##ELIF SMALL REDUCE
116 PARAMETER (IATBMX
= 8)
118 C..##IF LARGE XLARGE XXLARGE
120 PARAMETER (MAXHB
= 8000)
122 C..##ELIF REDUCE XSMALL
125 INTEGER MAXTRN
,MAXSYM
127 PARAMETER (MAXTRN
= 5000)
128 PARAMETER (MAXSYM
= 192)
131 C..##IF LONEPAIR (lonepair_max)
135 PARAMETER (MAXLP
= 2000)
136 PARAMETER (MAXLPH
= 4000)
138 C..##ENDIF (lonepair_max)
139 INTEGER NOEMAX
,NOEMX2
142 PARAMETER (NOEMAX
= 2000)
143 PARAMETER (NOEMX2
= 4000)
145 INTEGER MAXATC
, MAXCB
, MAXCH
, MAXCI
, MAXCP
, MAXCT
, MAXITC
, MAXNBF
148 PARAMETER (MAXATC
= 500, MAXCB
= 1500, MAXCH
= 3200, MAXCI
= 600,
149 & MAXCP
= 3000,MAXCT
= 15500,MAXITC
= 200, MAXNBF
=1000)
155 PARAMETER (MAXCN
= MAXITC*
(MAXITC
+1)/2)
156 INTEGER MAXA
, MAXAIM
, MAXB
, MAXT
, MAXP
157 INTEGER MAXIMP
, MAXNB
, MAXPAD
, MAXRES
158 INTEGER MAXSEG
, MAXGRP
159 C..##IF LARGE XLARGE XXLARGE
161 PARAMETER (MAXA
= SIZE
, MAXB
= SIZE
, MAXT
= SIZE
,
163 PARAMETER (MAXIMP
= 9200, MAXNB
= 17200, MAXPAD
= 8160,
167 PARAMETER (MAXSEG
= 1000)
176 PARAMETER (MAXAIM
= 2*SIZE
)
177 PARAMETER (MAXGRP
= 2*SIZE
/3)
179 INTEGER REDMAX
,REDMX2
182 PARAMETER (REDMAX
= 20)
183 PARAMETER (REDMX2
= 80)
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,
191 & MXRTX
= 5000, MXRTHA
= 300, MXRTHD
= 300,
193 & MXRTBL
= 5000, NICM
= 10)
194 INTEGER NMFTAB
, NMCTAB
, NMCATM
, NSPLIN
197 PARAMETER (NMFTAB
= 200, NMCTAB
= 3, NMCATM
= 12000, NSPLIN
= 3)
203 PARAMETER (MAXSHK
= SIZE*3
/4)
206 C..##IF IBM IBMRS CRAY INTEL IBMSP T3D REDUCE
208 PARAMETER (SCRMAX
= 5000)
214 PARAMETER (MXPIGG
=500)
216 INTEGER MXCOLO
,MXPUMB
217 PARAMETER (MXCOLO
=20,MXPUMB
=20)
220 INTEGER MAXUMP
, MAXEPA
, MAXNUM
223 PARAMETER (MAXUMP
= 10, MAXNUM
= 4)
227 PARAMETER (MAXING
=1000)
229 integer MAX_RINGSIZE
, MAX_EACH_SIZE
230 parameter (MAX_RINGSIZE
= 20, MAX_EACH_SIZE
= 1000)
232 parameter (MAXPATHS
= 8000)
233 integer MAX_TO_SEARCH
234 parameter (MAX_TO_SEARCH
= 6)
236 C-----------------------------------------------------------------------
237 C-----------------------------------------------------------------------
238 C:::##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
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
)
251 REAL*8 FIFTY
, SIXTY
, SVNTY2
, EIGHTY
, NINETY
, HUNDRD
,
252 & ONE2TY
, ONE8TY
, THRHUN
, THR6TY
, NINE99
, FIFHUN
, THOSND
,
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
)
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
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
)
282 PARAMETER (ANUM
=9999.0D0
, FMARK
=-999.0D0
)
283 PARAMETER (RSMALL
=1.0D
-10,RBIG
=1.0D20
)
289 C..##ELIF ALPHA T3D T3E
293 PARAMETER (RPRECI
= 2.22045D
-16, RBIGST
= 4.49423D
+307)
296 C-----------------------------------------------------------------------
297 C-----------------------------------------------------------------------
298 C:::##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
)
304 PARAMETER (COSMAX
=0.9999999999D0
)
306 PARAMETER (TIMFAC
=4.88882129D
-02)
308 PARAMETER (KBOLTZ
=1.987191D
-03)
313 PARAMETER (CCELEC
=332.0716D0
)
316 PARAMETER (CNVFRQ
=2045.5D0
/(2.99793D0*6
.28319D0
))
318 PARAMETER (SPEEDL
=2.99793D
-02)
320 PARAMETER (ATMOSP
=1.4584007D
-05)
322 PARAMETER (PATMOS
= 1.D0
/ ATMOSP
)
324 PARAMETER (BOHRR
= 0.529177249D0
)
326 PARAMETER (TOKCAL
= 627.5095D0
)
329 parameter(MDAKCAL
=143.9325D0
)
332 PARAMETER ( DEBYEC
= 2.541766D0
/ BOHRR
)
334 PARAMETER ( ZEROC
= 298.15D0
)
335 C-----------------------------------------------------------------------
336 C-----------------------------------------------------------------------
337 C:::##INCLUDE '~/charmm_fcm/exfunc.fcm'
342 CHARACTER*4 GTRMA
, NEXTA4
, CURRA4
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
,
352 * SRCHWD
, SRCHWS
, STRLNG
, DSIZE
, SSIZE
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
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
,
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
394 CHARACTER*8 ElementName
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
405 real*8 vangle
, OOPNGL
, TORNGL
, ElementMass
406 external vangle
, OOPNGL
, TORNGL
, ElementMass
408 C-----------------------------------------------------------------------
409 C-----------------------------------------------------------------------
410 C:::##INCLUDE '~/charmm_fcm/stack.fcm'
413 C...##IF LARGE XLARGE
414 C...##ELIF MEDIUM REDUCE
415 PARAMETER (STKSIZ
=4000000)
421 INTEGER LSTUSD
,MAXUSD
,STACK
422 COMMON /ISTACK
/ LSTUSD
,MAXUSD
,STACK
(STKSIZ
)
427 C-----------------------------------------------------------------------
428 C-----------------------------------------------------------------------
429 C:::##INCLUDE '~/charmm_fcm/heap.fcm'
431 C..##IFN UNICOS (unicos)
432 C...##IF XXLARGE (size)
433 C...##ELIF LARGE XLARGE (size)
434 C...##ELIF MEDIUM (size)
436 C....##ELIF TERRA (t3d2)
437 C....##ELIF ALPHA (t3d2)
438 C....##ELIF T3E (t3d2)
440 PARAMETER (HEAPDM
=2048000)
442 C...##ELIF SMALL (size)
443 C...##ELIF REDUCE (size)
444 C...##ELIF XSMALL (size)
447 INTEGER FREEHP
,HEAPSZ
,HEAP
448 COMMON /HEAPST
/ FREEHP
,HEAPSZ
,HEAP
(HEAPDM
)
449 LOGICAL LHEAP
(HEAPDM
)
450 EQUIVALENCE
(LHEAP
,HEAP
)
453 C..##IF SAVEFCM (save)
455 C-----------------------------------------------------------------------
456 C-----------------------------------------------------------------------
457 C:::##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
467 C-----------------------------------------------------------------------
468 C-----------------------------------------------------------------------
469 C:::##INCLUDE '~/charmm_fcm/deriv.fcm'
471 COMMON /DERIVR
/ DX
(MAXAIM
),DY
(MAXAIM
),DZ
(MAXAIM
)
474 C-----------------------------------------------------------------------
475 C-----------------------------------------------------------------------
476 C:::##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
,
484 & XTLTE
, XTLKE
, XTLPE
, XTLTEM
, XTLPEP
, XTLKEP
, XTLKP2
,
485 & TOT4
, TOTK4
, EPOT4
, TEM4
, MbMom
, BodyT
, PartT
487 & , SELF
, SCREEN
, COUL
,SOLV
, INTER
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,
501 & TOT4
= 37, TOTK4
= 38, EPOT4
= 39, TEM4
= 40,
502 & MbMom
= 41, BodyT
= 42, PartT
= 43
504 & , SELF
= 45, SCREEN
= 46, COUL
= 47,
505 & SOLV
= 48, INTER
= 49
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
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,
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,
569 & VIXX
= 10, VIXY
= 11, VIXZ
= 12, VIYX
= 13,
570 & VIYY
= 14, VIYZ
= 15, VIZX
= 16, VIZY
= 17,
572 & PEXX
= 19, PEXY
= 20, PEXZ
= 21, PEYX
= 22,
573 & PEYY
= 23, PEYZ
= 24, PEZX
= 25, PEZY
= 26,
575 & PIXX
= 28, PIXY
= 29, PIXZ
= 30, PIYX
= 31,
576 & PIYY
= 32, PIYZ
= 33, PIZX
= 34, PIZY
= 35,
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
)
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
)
595 INTEGER ECALLS
, TOT1ST
, TOT2ND
596 COMMON /EMISCI
/ ECALLS
, TOT1ST
, TOT2ND
597 REAL*8 EOLD
, FITA
, DRIFTA
, EAT0A
, CORRA
, FITP
, DRIFTP
,
599 COMMON /EMISCR
/ EOLD
, FITA
, DRIFTA
, EAT0A
, CORRA
,
600 & FITP
, DRIFTP
, EAT0P
, CORRP
614 REAL*8 TSMTRM
(LENENT
),TSMTMP
(LENENT
)
615 COMMON /TSMENG
/ TSMTRM
,TSMTMP
624 C-----------------------------------------------------------------------
625 C-----------------------------------------------------------------------
626 C:::##INCLUDE '~/charmm_fcm/dimb.fcm'
627 C..##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
649 C-----------------------------------------------------------------------
650 C-----------------------------------------------------------------------
651 C:::##INCLUDE '~/charmm_fcm/ctitla.fcm'
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
)
660 C-----------------------------------------------------------------------
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
(*)
674 LOGICAL LNOMA
,LRAISE
,LSCI
,LBIG
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
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
690 LOGICAL LCARD
,LAPPE
,LPURG
,LWDINI
,QCALC
,QMASWT
,QMIX
,QDIAG
706 NFREG6
=(NFREG
-6)/NPAR
709 IF(NFREG
.GT
.PARDIM
) CALL WRNDIE
(-3,'<NMDIMB>',
710 1 'NFREG IS LARGER THAN PARDIM*3')
712 C ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
716 C ALLOCATE-SPACE-FOR-DIAGONALIZATION
720 C ALLOCATE-SPACE-FOR-REDUCED-BASIS
724 C ALLOCATE-SPACE-FOR-OTHER-ARRAYS
729 C Space allocation for working arrays of EISPACK
730 C diagonalization subroutines
732 C ALLOCATE-SPACE-FOR-LSCI
737 C ALLOCATE-DUMMY-SPACE-FOR-LSCI
744 LENCM
=INBCMP
(NATOM
-1)*9+NATOM*6
751 CALL ENERGY
(X
,Y
,Z
,DX
,DY
,DZ
,BNBND
,BIMAG
,NAT3
,DD1CMP
,.TRUE
.,1)
755 C Mass weight DD1CMP matrix
757 CALL MASSDD
(DD1CMP
,DDM
,INBCMP
,JNBCMP
,NATOM
)
759 CALL WRNDIE
(-3,'<NMDIMB>','QDISK OPTION NOT SUPPORTED YET')
767 C Fill DDV with six translation-rotation vectors
769 CALL TRROT
(X
,Y
,Z
,DDV
,NAT3
,1,DDM
)
770 CALL CPARAY
(HEAP
(TRAROT
),DDV
,NAT3
,1,6,1)
774 CALL ORTHNM
(1,6,NTR
,HEAP
(TRAROT
),NAT3
,.FALSE
.,TOLER
)
776 IF(IUNRMD
.LT
. 0) THEN
778 C If no previous basis is read
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
/)
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
),
794 IF(PRNLEV
.GE
.9) CALL PRINTE
(OUTU
,EPROP
,ETERM
,'VIBR',
795 1 'ENR',.TRUE
.,1,ZERO
,ZERO
)
797 C Generate the lower section of the matrix and diagonalize
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
)
815 C Put the PARDDV vectors into DDV and replace the elements which do
816 C not belong to the considered partitioned region by zeros.
818 CALL ADJNME
(DDV
,PARDDV
,NAT3
,NDIM
,NFRE
,NFRET
,IS1
,IS2
)
821 PARDDF
(J
)=CNVFRQ*SQRT
(ABS
(PARDDE
(J
)))
822 IF(PARDDE
(J
) .LT
. 0.0) PARDDF
(J
)=-PARDDF
(J
)
827 PARDDF
(J
)=CNVFRQ*SQRT
(ABS
(PARDDE
(J
)))
828 IF(PARDDE
(J
) .LT
. 0.0) PARDDF
(J
)=-PARDDF
(J
)
834 WRITE(OUTU
,516) (J
,PARDDF
(J
),J
=1,NFRE
)
837 IF(NFRET
.GE
. NFREG
) GOTO 10
839 512 FORMAT(/' NMDIMB: Diagonalization of part',I5
,' completed')
840 514 FORMAT(' NMDIMB: Frequencies'/)
841 516 FORMAT(5(I4
,F12
.6
))
844 C Orthonormalize the eigenvectors
848 CALL ORTHNM
(1,NFRET
,NFRET
,DDV
,NAT3
,LPURG
,TOLER
)
851 C Do reduced basis diagonalization using the DDV vectors
852 C and get eigenvectors of zero iteration
856 WRITE(OUTU
,523) NFRET
858 521 FORMAT(/' NMDIMB: Iteration number = ',I5
)
859 523 FORMAT(' NMDIMB: Dimension of the reduced basis set = ',I5
)
861 IF(PRNLEV
.GE
.2) WRITE(OUTU
,585) NFRET
,IUNMOD
862 525 FORMAT(' NMDIMB: ',I5
,' basis vectors are saved in unit',I5
)
865 CALL WRTNMD
(LCARD
,1,NFRET
,NAT3
,DDV
,DDSCR
,DDEV
,IUNMOD
,AMASS
)
868 CALL CPARAY
(HEAP
(DDVBAS
),DDV
,NAT3
,1,NFRET
,1)
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
)
878 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
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
)
897 C Read in existing basis
901 531 FORMAT(/' NMDIMB: Calculations restarted')
908 CALL RDNMD
(LCARD
,NFRET
,NFREG
,NAT3
,NDIM
,
909 1 DDV
,DDSCR
,DDF
,DDEV
,
910 2 IUNRMD
,LAPPE
,ISTRT
,ISTOP
)
912 IF(NFRET
.GT
.NFREG
) THEN
914 CALL WRNDIE
(-1,'<NMDIMB>',
915 1 'Not enough space to hold the basis. Increase NMODes')
919 WRITE(OUTU
,533) NFRET
,IUNRMD
921 WRITE(OUTU
,516) (J
,DDF
(J
),J
=1,NFRET
)
923 533 FORMAT(/' NMDIMB: ',I5
,' restart modes read from unit ',I5
)
927 C -------------------------------------------------
928 C Here starts the mixed-basis diagonalization part.
929 C -------------------------------------------------
932 C Check cut-off frequency
934 CALL SELNMD
(DDF
,NFRET
,CUTF1
,NFCUT1
)
937 IF(NFCUT1*2
-6.GT
.NFREG
) THEN
938 IF(PRNLEV
.GE
.2) WRITE(OUTU
,537) DDF
(NFRRES
)
945 537 FORMAT(/' NMDIMB: Too many vectors for the given cutoff frequency'
946 1 /' Cutoff frequency is decreased to',F9
.3
)
948 C Compute the new partioning of the molecule
950 CALL PARTIC
(NAT3
,NFREG
,NFCUT1
,NPARMX
,NPARC
,ATMPAR
,NFRRES
,
954 ATMPAS
(1,I
)=ATMPAR
(1,I
)
955 ATMPAS
(2,I
)=ATMPAR
(2,I
)
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
.
966 WRITE(OUTU
,543) ITER
,ITMX
967 IF(QDW
) WRITE(OUTU
,545) IPAR1
,IPAR2
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
)
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
,
978 C If double windowing is defined, the original block sizes are divided
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
)
994 C Find out with which block to continue (double window method only)
1005 CALL IPART
(HEAP
(SUBLIS
),II
,IPAR1
,IPAR2
,HEAP
(ATMPAF
),
1007 IF((IPAR1
.EQ
.IPA1
).AND
.(IPAR2
.EQ
.IPA2
)) GOTO 500
1014 DO WHILE((CVGMX
.GT
.TOLDIM
).AND
.(ITER
.LT
.ITMX
))
1017 IF(PRNLEV
.GE
.2) WRITE(OUTU
,553) ITER
1018 553 FORMAT(/' NMDIMB: Iteration number = ',I8
)
1024 CALL PARTDS
(NAT3
,NPARC
,ATMPAR
,NPARS
,ATMPAS
,INIDS
,NPARMX
,
1025 1 DDF
,NFREG
,CUTF1
,PARDIM
,NFCUT1
)
1026 C DO-THE-DIAGONALISATIONS
1031 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1043 CALL IPART
(HEAP
(SUBLIS
),II
,IPAR1
,IPAR2
,HEAP
(ATMPAF
),
1048 IF(PRNLEV
.GE
.2) WRITE(OUTU
,553) ITER
1049 C DO-THE-DWIN-DIAGONALISATIONS
1054 IF((IRESF
.EQ
.SAVF
).OR
.(ITER
.EQ
.ITMX
)) THEN
1057 C DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1062 IF((CVGMX
.LE
.TOLDIM
).OR
.(ITER
.EQ
.ITMX
)) GOTO 600
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
)
1083 C-----------------------------------------------------------------------
1084 C INTERNAL PROCEDURES
1085 C-----------------------------------------------------------------------
1086 C TO DO-THE-DIAGONALISATIONS-WITH-RESIDUALS
1088 IF(IUNRMD
.LT
.0) THEN
1089 CALL SELNMD
(DDF
,NFRET
,CUTF1
,NFC
)
1093 IF(NFCUT*2
-6 .GT
. NFREG
) THEN
1096 IF(PRNLEV
.GE
.2) THEN
1097 WRITE(OUTU
,562) ITER
1098 WRITE(OUTU
,564) CUTF1
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
)
1109 C do reduced diagonalization with preceding eigenvectors plus
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
)
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
)
1127 CALL ORTHNM
(1,NFRET
,NFCUT
,DDV
,NAT3
,LPURG
,TOLER
)
1130 IF(PRNLEV
.GE
.2) WRITE(OUTU
,568) NFRET
1131 568 FORMAT(' after orthonormalization = ',I5
)
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
)
1137 CALL WRTNMD
(LCARD
,1,NFCUT
,NAT3
,DDV
,DDSCR
,DDEV
,IUNMOD
,AMASS
)
1140 CALL CPARAY
(HEAP
(DDVBAS
),DDV
,NAT3
,1,NFCUT
,1)
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
)
1154 C-----------------------------------------------------------------------
1155 C TO DO-THE-DIAGONALISATIONS
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')
1168 IF(NFRET
.GT
.NFREG
) NFRET
=NFREG
1169 CALL CLETR
(DDV
,HEAP
(TRAROT
),NAT3
,1,NFCUT1
,NFCUT
,DDEV
,DDF
)
1171 CALL ADZER
(DDV
,1,NFCUT1
,NAT3
,IS1
,IS2
)
1175 CALL ORTHNM
(1,NFCUT1
,NFCUT
,DDV
,NAT3
,LPURG
,TOLER
)
1177 CALL CPARAY
(HEAP
(DDVBAS
),DDV
,NAT3
,1,NFCUT
,1)
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
)
1188 IF(NFCUT
.GT
.NFRRES
) NFCUT
=NFRRES
1194 C-----------------------------------------------------------------------
1195 C TO DO-THE-DWIN-DIAGONALISATIONS
1198 C Store the DDV vectors into DDVBAS
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 ',
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')
1213 IF(NFRET
.GT
.NFREG
) NFRET
=NFREG
1215 C Prepare the DDV vectors consisting of 6 translations-rotations
1216 C + eigenvectors from 7 to NFCUT1 + cartesian displacements vectors
1217 C spanning the atoms from IS1 to IS2
1219 CALL CLETR
(DDV
,HEAP
(TRAROT
),NAT3
,1,NFCUT1
,NFCUT
,DDEV
,DDF
)
1222 CALL ADZERD
(DDV
,1,NFCUT1
,NAT3
,IS1
,IS2
,IS3
,IS4
)
1225 CALL ORTHNM
(1,NFCUT1
,NFCUT
,DDV
,NAT3
,LPURG
,TOLER
)
1227 CALL CPARAY
(HEAP
(DDVBAS
),DDV
,NAT3
,1,NFCUT
,1)
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
)
1240 IF(NFCUT
.GT
.NFRRES
) NFCUT
=NFRRES
1245 C-----------------------------------------------------------------------
1248 IF(PRNLEV
.GE
.2) WRITE(OUTU
,583) IUNMOD
1249 583 FORMAT(/' NMDIMB: Saving the eigenvalues and eigenvectors to unit'
1251 REWIND
(UNIT
=IUNMOD
)
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
,
1262 C-----------------------------------------------------------------------
1263 C TO ALLOCATE-SPACE-FOR-DIAGONALIZATION
1265 DDV2
=ALLHP
(IREAL8
((PARDIM
+3)*(PARDIM
+3)))
1266 JSPACE
=IREAL8
((PARDIM
+4))*8
1267 JSP
=IREAL8
(((PARDIM
+3)*(PARDIM
+4))/2)
1273 C-----------------------------------------------------------------------
1274 C TO ALLOCATE-SPACE-FOR-REDUCED-BASIS
1277 DDVBAS
=ALLHP
(IREAL8
(NAT3
))
1279 DDVBAS
=ALLHP
(IREAL8
(NFREG*NAT3
))
1283 C-----------------------------------------------------------------------
1284 C TO ALLOCATE-SPACE-FOR-TRANSROT-VECTORS
1286 TRAROT
=ALLHP
(IREAL8
(6*NAT3
))
1289 C-----------------------------------------------------------------------
1290 C TO ALLOCATE-SPACE-FOR-LSCI
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))
1304 C-----------------------------------------------------------------------
1305 C TO ALLOCATE-DUMMY-SPACE-FOR-LSCI
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))
1319 C-----------------------------------------------------------------------
1320 C TO ALLOCATE-SPACE-FOR-OTHER-ARRAYS
1322 IUPD
=ALLHP
(INTEG4
(PARDIM
+3))