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