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