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