1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck srestr */ 20 SUBROUTINE SRESTR 21C 22C Several SIRIUS-variables has to be restored, because they are overwritten 23C by ABACUS in SETSIR module. Called from SIRRDI. 24C 25#include "implicit.h" 26#include "maxorb.h" 27#include "gnrinf.h" 28#include "infinp.h" 29#include "inpbak.h" 30#include "infpri.h" 31#include "pribak.h" 32#include "pcmlog.h" 33C 34C Restoring should only be done if SIRIUS has been through its 35C input processing. 36C 37 IF (.NOT. SIR_INPPRC) RETURN 38C 39C /INFINP/ values are restored 40C 41 ISPIN = ISPINB 42 ISTATE = ISTATB 43 LSYM = LSYMB 44 NACTEL = NACTLB 45 MCTYPE = MCTYPB 46 LSOLMX = LSLMXB 47 NLMSOL = NLMSLB 48 NELMN1 = NLMN1B 49 NELMX1 = NLMX1B 50 NELMN3 = NLMN3B 51 NELMX3 = NLMX3B 52 LROOTS = LROTSB 53 NROOTS = NROTSB 54 IORTO = IORTOB 55 ICI0 = ICI0B 56 KDEL = KDELB 57 ICHECK = ICHCKB 58 NTIT = NTITB 59 MAXMAC = MXMACB 60 MAXMIC = MXMICB 61 MAXJT = MAXJTB 62 MAXCIT = MXCITB 63 MAXUIT = MXUITB 64 MAXAPM = MXAPMB 65 MAXABS = MXABSB 66 ITRLVL = ITRLVB 67 ITRFIN = ITRFNB 68 JCHSYM = JCHSMB 69 JCHORB = JCHRBB 70 NROOCI = NROCIB 71 ISTACI = ISTCIB 72 NFIELD = NFILDB 73 MXCIMA = MXCMAB 74 ICICNO = ICCNOB 75 IMCCNO = IMCNOB 76 IF (NEWSYM) THEN 77 DOSCF = DOSCFB 78 DOMP2 = DOMP2B 79 DOCINO = DOCNOB 80 DOCI = DOCIB 81 FCVORB = FCVRBB 82 ELSE 83 DOSCF = (DOSCFB .AND. .NOT. DOMCB) 84 DOMP2 = (DOMP2B .AND. .NOT. DOMCB) 85 DOCINO = (DOCNOB .AND. .NOT. DOMCB) 86 DOCI = (DOCIB .AND. .NOT. DOMCB) 87 FCVORB = (FCVRBB .AND. .NOT. DOMCB) 88 END IF 89 DOMC = DOMCB 90 DORSP = DORSPB 91 LNOROT = LNROTB 92 LMOORD = LMORDB 93 DIRFCK = DRFCKB 94 CORHOL = CORHLB 95 CORRLX = CRRLXB 96 RESPHP = RSPHPB 97 JOLSEN = JOLSNB 98 ABAIPH = ABIPHB 99 INERSI = INRSIB 100 INERSF = INRSFB 101 SUPSYM = SUPSMB 102 SPIN = SPINB 103 POTNUC = POTNCB 104 EPSOL = EPSOLB 105 EPSTAT = EPSTTB 106 PCM = PCMB 107 RSOL(1)= RSOLB(1) 108 RSOL(2)= RSOLB(2) 109 RSOL(3)= RSOLB(3) 110 THRGRD = TRGRDB 111 THRPWF = TRPWFB 112 THRCI = TRCIB 113 THRMC = TRMCB 114 THRCGR = TRCGRB 115 CMAXMO = CMXMOB 116 THROVL = TROVLB 117 THRSSY = TRSSYB 118 DO I = 1, MAXRTS 119 IROOT(I) = IROOTB(I) 120 END DO 121 DO I = 1, MXCORB 122 NOROT(I) = NOROTB(I) 123 IMOORD(I) = IMORDB(I) 124 CENT(I) = CENTB(I) 125 TYPE(I) = TYPEB(I) 126 END DO 127 DO I = 1, NFLAG 128 FLAG(I) = FLAGB(I) 129 END DO 130 DO I = 1, 6 131 TITLE(I) = TITLEB(I) 132 END DO 133 TITMOL(1) = TITMLB(1) 134 TITMOL(2) = TITMLB(2) 135 DO I = 1, MXFELT 136 EFIELD(I) = EFILDB(I) 137 LFIELD(I) = LFILDB(I) 138 END DO 139C 140C /INFPRI/ values are restored 141C 142 DO I = 1, NPFLAG 143 P4FLAG(I) = P4FLGB(I) 144 P6FLAG(I) = P6FLGB(I) 145 END DO 146 IPRI4 = IPRI4B - 2 147 IPRI6 = IPRI6B 148 IPRSIR = IPRSRB 149 IPRCNO = IPRCNB 150 IPRDIA = IPRDIB 151 IPRSIG = IPRSGB 152 IPRDNS = IPDNSB 153 IPRSOL = IPRSLB 154 IPRKAP = IPRKPB 155 IPRMUL = IPRMLB 156 IPRCIX = IPRCXB 157 IPRRHF = IPRHFB 158 IPRFCK = IPRFCB 159 IPRAVE = IPRAVB 160 MPRI4 = MPRI4B 161 MPRI6 = MPRI6B 162 MPRSIR = MPRSRB 163 RETURN 164 END 165C /* Deck ssave */ 166 SUBROUTINE SSAVE 167C 168C Several SIRIUS-variables has to be saved, because they are overwritten 169C by ABACUS in SETSIR module. Called from SIRRDI. 170C 171#include "implicit.h" 172#include "maxorb.h" 173#include "gnrinf.h" 174#include "infinp.h" 175#include "inpbak.h" 176#include "infpri.h" 177#include "pribak.h" 178#include "pcmlog.h" 179C 180C /INFINP/ values are saved to /INFBAK/ 181C 182 ISPINB = ISPIN 183 ISTATB = ISTATE 184 LSYMB = LSYM 185 NACTLB = NACTEL 186 MCTYPB = MCTYPE 187 LSLMXB = LSOLMX 188 NLMSLB = NLMSOL 189 NLMN1B = NELMN1 190 NLMX1B = NELMX1 191 NLMN3B = NELMN3 192 NLMX3B = NELMX3 193 LROTSB = LROOTS 194 NROTSB = NROOTS 195 IORTOB = IORTO 196 ICI0B = ICI0 197 KDELB = KDEL 198 ICHCKB = ICHECK 199 NTITB = NTIT 200 MXMACB = MAXMAC 201 MXMICB = MAXMIC 202 MAXJTB = MAXJT 203 MXCITB = MAXCIT 204 MXUITB = MAXUIT 205 MXAPMB = MAXAPM 206 MXABSB = MAXABS 207 ITRLVB = ITRLVL 208 ITRFNB = ITRFIN 209 JCHSMB = JCHSYM 210 JCHRBB = JCHORB 211 NROCIB = NROOCI 212 ISTCIB = ISTACI 213 NFILDB = NFIELD 214 MXCMAB = MXCIMA 215 ICCNOB = ICICNO 216 IMCNOB = IMCCNO 217 DOSCFB = DOSCF 218 DOMP2B = DOMP2 219 DOCNOB = DOCINO 220 DOCIB = DOCI 221 DOMCB = DOMC 222 DORSPB = DORSP 223 FCVRBB = FCVORB 224 LNROTB = LNOROT 225 LMORDB = LMOORD 226 DRFCKB = DIRFCK 227 CORHLB = CORHOL 228 CRRLXB = CORRLX 229 RSPHPB = RESPHP 230 JOLSNB = JOLSEN 231 ABIPHB = ABAIPH 232 INRSIB = INERSI 233 INRSFB = INERSF 234 SUPSMB = SUPSYM 235 SPINB = SPIN 236 POTNCB = POTNUC 237 EPSOLB = EPSOL 238 EPSTTB = EPSTAT 239 PCMB = PCM 240 RSOLB(1)= RSOL(1) 241 RSOLB(2)= RSOL(2) 242 RSOLB(3)= RSOL(3) 243 TRGRDB = THRGRD 244 TRPWFB = THRPWF 245 TRCIB = THRCI 246 TRMCB = THRMC 247 TRCGRB = THRCGR 248 CMXMOB = CMAXMO 249 TROVLB = THROVL 250 TRSSYB = THRSSY 251 TITLEB(1) = TITLE(1) 252 TITLEB(2) = TITLE(2) 253 TITLEB(3) = TITLE(3) 254 TITLEB(4) = TITLE(4) 255 TITLEB(5) = TITLE(5) 256 TITLEB(6) = TITLE(6) 257 DO 100 I = 1, MAXRTS 258 IROOTB(I) = IROOT(I) 259 100 CONTINUE 260 DO 110 I = 1, MXCORB 261 NOROTB(I) = NOROT(I) 262 IMORDB(I) = IMOORD(I) 263 CENTB(I) = CENT(I) 264 TYPEB(I) = TYPE(I) 265 110 CONTINUE 266 DO 130 I = 1, NFLAG 267 FLAGB(I) = FLAG(I) 268 130 CONTINUE 269 TITMLB(1) = TITMOL(1) 270 TITMLB(2) = TITMOL(2) 271 DO 150 I = 1, MXFELT 272 EFILDB(I) = EFIELD(I) 273 LFILDB(I) = LFIELD(I) 274 150 CONTINUE 275C 276C /INFPRI/ values are saved 277C 278 DO I = 1, NPFLAG 279 P4FLGB(I) = P4FLAG(I) 280 P6FLGB(I) = P6FLAG(I) 281 END DO 282 IPRI4B = IPRI4 283 IPRI6B = IPRI6 284 IPRSRB = IPRSIR 285 IPRCNB = IPRCNO 286 IPRDIB = IPRDIA 287 IPRSGB = IPRSIG 288 IPDNSB = IPRDNS 289 IPRSLB = IPRSOL 290 IPRKPB = IPRKAP 291 IPRMLB = IPRMUL 292 IPRCXB = IPRCIX 293 IPRHFB = IPRRHF 294 IPRFCB = IPRFCK 295 IPRAVB = IPRAVE 296 MPRI4B = MPRI4 297 MPRI6B = MPRI6 298 MPRSRB = MPRSIR 299 RETURN 300 END 301C /* Deck prestr */ 302 SUBROUTINE PRESTR 303C 304C ABACUS-variables concerning execution of modules are restored. 305C Called from ABADRV. Necessary because modules should be run 306C in each geometry iteration. 307C 308#include "implicit.h" 309#include "past.h" 310#include "pastbk.h" 311 PASTWO = PSTWOB 312 PASORT = PSORTB 313 PASDIP = PSDIPB 314 PASONE = PSONEB 315 PASRES = PSRESB 316 PASREL = PSRELB 317 PASDPR = PSDPRB 318 PASCRS = PSCRSB 319 PASCZR = PSCZRB 320 PASCTR = PSCTRB 321 PASCRL = PSCRLB 322 PASMAG = PSMAGB 323 PASAAT = PSAATB 324 PASRTR = PSRTRB 325 PASLRS = PSLRSB 326 PASTRP = PSTRPB 327 PASLNR = PSLNRB 328 PASEXC = PSEXCB 329 RETURN 330 END 331C /* Deck psave */ 332 SUBROUTINE PSAVE 333C 334C ABACUS-variables concerning execution of modules are saved. 335C Called from ABADRV. 336C 337#include "implicit.h" 338#include "past.h" 339#include "pastbk.h" 340 PSTWOB = PASTWO 341 PSORTB = PASORT 342 PSDIPB = PASDIP 343 PSONEB = PASONE 344 PSRESB = PASRES 345 PSRELB = PASREL 346 PSDPRB = PASDPR 347 PSCRSB = PASCRS 348 PSCZRB = PASCZR 349 PSCTRB = PASCTR 350 PSCRLB = PASCRL 351 PSMAGB = PASMAG 352 PSAATB = PASAAT 353 PSRTRB = PASRTR 354 PSLRSB = PASLRS 355 PSTRPB = PASTRP 356 PSLNRB = PASLNR 357 PSEXCB = PASEXC 358 RETURN 359 END 360C /* Deck bndchk */ 361 SUBROUTINE BNDCHK(WORK, LMWORK2, WRKDLM, PROG) 362#include "implicit.h" 363 DIMENSION WORK(LMWORK2) 364 CHARACTER*(*) PROG 365#include "priunit.h" 366C 367C Check memory traps. Gives error message if any of the programs 368C have been outside the declared memory area. 369C 370 IF (WORK(1) .NE. WRKDLM .OR. 371 & WORK(LMWORK2) .NE. WRKDLM) THEN 372 WRITE (LUPRI,'(//3A)') 373 * ' WARNING, ',PROG,' module has been out of bounds.' 374 IF (WORK(1) .NE. WRKDLM) WRITE (LUPRI,'(/A)') 375 * ' WORK(0) has been destroyed' 376 IF (WORK(LMWORK2) .NE. WRKDLM) WRITE (LUPRI,'(/A)') 377 * ' WORK(LMWORK+1) has been destroyed' 378 CALL QUIT('WARNING, ' // PROG // ' has been out of bounds.') 379 END IF 380 RETURN 381 END 382C /* Deck exeher */ 383 SUBROUTINE EXEHER(WORK, LMWORK, WRKDLM) 384 385 use qmcmm, only: read_pot_qmnpmm 386 387#include "implicit.h" 388#include "dummy.h" 389#include "maxaqn.h" 390#include "mxcent.h" 391 DIMENSION WORK(0:LMWORK) 392#include "priunit.h" 393#include "gnrinf.h" 394! cbiher.h : TWOBYTEPACKING 395! nuclei.h : NBASIS 396#include "cbiher.h" 397#include "nuclei.h" 398#include "dftcom.h" 399#include "pcmlog.h" 400#include "qm3.h" 401C 402C Run Integral section 403C 404 CALL QENTER('HERMIT') 405 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini HERMIT') 406 NEWPRP = .TRUE. 407 NEWMAT = .TRUE. 408 DFTGRID_DONE = .FALSE. 409 DFTGRID_DONE_OLD = .FALSE. 410C ... we need new grid for numerical integrations when new geometry 411C 412C *************************************************************** 413C Reading in the interaction parameters for the MM system types: 414C *************************************************************** 415C 416 IF (QM3 .AND. .NOT.ONLYOV) CALL QM3_RDPOT() 417C 418C *************************************************************** 419C Reading in the interaction parameters for the NP and MM 420C subsystems needed by QM/NP/MM embedding scheme 421C *************************************************************** 422C 423 IF (QMNPMM .AND. .NOT.ONLYOV) THEN 424 CALL READ_POT_QMNPMM() 425 END IF 426C 427C *************************************************************** 428C Read initial geometry, orbital spec., etc. (MOLECULE format). 429C *************************************************************** 430C 431C If we use the Douglas-Kroll transformation, we first generate 432C the necessary integrals in the primitive basis. 433C See further notes in DKH_INTF routine. 434C 435 IF (DKTRAN .AND. RNHERM) THEN 436 DKHINT = .TRUE. ! now we are calculating DKH 1-el. integrals 437 CALL DKH_INTF(WORK(1),LMWORK) 438 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'DKH_INTF ') 439 DKHINT = .FALSE. 440 NEWPRP = .TRUE. 441 END IF 442C 443C Third parameter .TRUE. in call indicates that LUONEL will be written. 444C 445 CALL READIN(WORK(1),LMWORK,RNHERM) 446 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'READIN') 447C 448C Initialize Hermit variables, also when RNHERM false (e.g. for ABACUS) 449C 450 CALL HERINI 451 TWOBYTEPACKING = TWOBYTEPACKING .OR. NBASIS .GT. 255 452 ! we need two bytes to pack indices if NBASIS .gt. 255; we also need 453 ! the .OR. test, because TWOBYTEPACKING might have been requested in input 454C 455C *************************************************************** 456C Setup of MM information. 457C *************************************************************** 458C 459 IF (QM3 .AND. .NOT.ONLYOV) CALL QM3_SETUP() 460 CALL FLSHFO(LUPRI) 461C 462C *************************************************************** 463C Now we generate all requested integrals in the contracted basis 464C (if DKTRAN: including transformation of the Douglas-Kroll 465C integrals to contracted basis). 466C 467 IF (RNHERM) THEN 468 CALL TITLER('Starting in Integral Section (HERMIT)',' ',200) 469 CALL HERCTL(WORK(1),LMWORK) 470C ... Note: if RUNTWO or RUNERI then HERCTL sets FTRCTL = .TRUE. 471C for transformation control, old AO/MO files cannot be used any more. 472C *************************************************************** 473 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'HERMIT') 474 CALL TITLER('End of Integral Section (HERMIT)',' ',200) 475 CALL FLSHFO(LUPRI) 476 END IF 477 CALL QEXIT('HERMIT') 478 RETURN 479 END 480C /* Deck exesir */ 481 SUBROUTINE EXESIR(WORK, LMWORK, WRKDLM) 482#include "implicit.h" 483 DIMENSION WORK(0:LMWORK) 484#include "priunit.h" 485#include "maxorb.h" 486#include "gnrinf.h" 487#include "inftra.h" 488#include "cbihr2.h" 489#include "inftap.h" 490#include "huckel.h" 491C 492C 493C Run SIRIUS 494C 495 CALL QENTER('SIRIUS') 496C 497 CALL TITLER('Starting in Wave Function Section (SIRIUS)',' ',200) 498 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini SIRIUS') ! boundary checks on WORK 499C 500C SIRINI defines buffer lengths for I/O and the IROW() array 501C 502 CALL SIRINI ! SIRius INItialization 503 WRINDX = .TRUE. ! new sirius call, make sure orbital rotation info is written to file in SETWOP 504C 505 CALL SIRCTL(ICONV,WORK(1),LMWORK) 506 ! ICONV > 0: converged 507 ! = 0: not converged 508 ! < 0: convergence not requested (e.g. restart from old MCSCF) 509 510 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'SIRIUS') ! boundary checks on WORK 511C 512C If not numerical differentiation: 513C We should now have MO's on file, no need for DOHUCKEL anymore 514C 515CHJ IF (.NOT.NMWALK) THEN 516CHJ: May09: we try to set DOHUCKEL true later if NEWSYM true instead 517CHJ HUCKEL guess is only used in SIRIUS if first call or NEWSYM 518 DOHUCKEL = .FALSE. 519CHJ END IF 520C 521C Several files left open, have to be closed manually 522C 523C With the present version of the code, AOSUPINT cannot be used either 524C in RESPONSE or ABACUS, and for this reason we delete the file here 525C 526 IF (LUSUPM.LE.0) CALL GPOPEN(LUSUPM, 527 & FNSUPM,'UNKNOWN',' ','UNFORMATTED',IDUMMY,.FALSE.) 528 CALL GPCLOSE(LUSUPM,'DELETE') 529 IF (LUIT2 .GT. 0) CALL GPCLOSE(LUIT2,'DELETE') 530 IF (LUIT3 .GT. 0) CALL GPCLOSE(LUIT3,'DELETE') 531 IF (LUIT5 .GT. 0) CALL GPCLOSE(LUIT5,'DELETE') 532 IF (LUINTM .GT. 0) CALL GPCLOSE(LUINTM,'KEEP') 533 CALL TITLER('End of Wave Function Section (SIRIUS)',' ',200) 534 535 IF (ICONV .EQ. 0 .AND. .NOT.WLKREJ) THEN 536 WRITE(LUPRI, '(/A)') 537 & '@ - DALTON aborted because wave function not converged!' 538 CALL QUIT('DALTON aborted because wave function not converged') 539 END IF 540 CALL FLSHFO(LUPRI) 541 CALL QEXIT('SIRIUS') 542 RETURN 543 END 544C /* Deck execc */ 545 SUBROUTINE EXECC(WORK, LMWORK, WRKDLM) 546#include "implicit.h" 547 DIMENSION WORK(0:LMWORK) 548#include "priunit.h" 549#include "gnrinf.h" 550#include "maxorb.h" 551#include "mxcent.h" 552#include "inftra.h" 553#include "optinf.h" 554#include "nuclei.h" 555C 556C Run Coupled Cluster program 557C 558 CALL QENTER('CC') 559C 560 CALL TITLER('Starting in Coupled Cluster Section (CC)',' ',200) 561 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini CC') 562C 563 IF (RNSIRI .AND. .NOT. RNHERM) THEN 564 CALL SETHER(IPRUSR,.FALSE.,WORK(1),LMWORK) 565 CALL ER2INI 566 END IF 567 568! 569! The CC module uses a lot of static allocation for pointer arrays. 570! To minimize the static allocation in Dalton we use a smaller 571! max number of orbitals in CC than the global MXCORB. /hjaaj Sep 2013 572! 573 IF (NBASIS .GT. MXCORB_CC) THEN 574 WRITE (LUPRI,'(//A,I6/A,I6/A)') 575 & ' ERROR Too many orbitals for CC module:',NBASIS, 576 & ' ERROR Limit is MXCORB_CC in maxorb.h :',MXCORB_CC, 577 % ' ERROR To proceed increase MXCORB_CC in maxorb.h and recompile' 578 CALL QUIT('Too many orbitals for CC module.'// 579 & 'Increase MXCORB_CC in maxorb.h and recompile.') 580 END IF 581 582 ! ITRNMR > 0 required for proper restart of CPHF solver 583 ISAVITR = ITRNMR 584 ITRNMR = 1 585 586 CALL CC_DRV(WORK(1),LMWORK) 587 588 ITRNMR = ISAVITR 589 590 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'CC') 591 592 CALL TITLER('End of Coupled Cluster Section (CC)',' ',200) 593 CALL FLSHFO(LUPRI) 594 CALL QEXIT('CC') 595 RETURN 596 END 597C /* Deck exeaba */ 598 SUBROUTINE EXEABA(WORK, LMWORK, WRKDLM) 599 600 use pelib_interface, only: use_pelib, pelib_ifc_activate, 601 & pelib_ifc_deactivate 602 603#include "implicit.h" 604 DIMENSION WORK(0:LMWORK) 605#include "priunit.h" 606#include "maxorb.h" 607#include "mxcent.h" 608#include "inftra.h" 609#include "inftap.h" 610#include "huckel.h" 611#include "gnrinf.h" 612#include "abainf.h" 613#include "pcmlog.h" 614 615 LOGICAL SAVE_PELIB 616C 617C Run ABACUS 618C 619 IF (SKIPAB) RETURN 620 CALL QENTER('ABACUS') 621 CALL TITLER('Starting in Static Property Section (ABACUS) -', 622 & ' ',200) 623 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini ABACUS') 624 625C Test if a requested molecular derivative calculation is implemented 626C NAORDR = max order of analytical geometry derivatives 627 628 SAVE_PELIB = .FALSE. 629 IF (USE_PELIB()) THEN 630 SAVE_PELIB = .TRUE. 631 CALL PELIB_IFC_DEACTIVATE() 632 END IF 633 CALL FNDANA(NAORDR) 634 IF (SAVE_PELIB) THEN 635 IF (MOLHES .AND. NAORDR .EQ. 2) THEN 636 WRITE(LUPRI,'(//A)') 637 & '@ WARNING. Analytical molecular Hessian is approximate'// 638 & ' because embedding contributions are neglected.' 639 END IF 640 CALL PELIB_IFC_ACTIVATE() 641 END IF 642 IF (MOLHES .AND. NAORDR .LT. 2) THEN 643 WRITE(LUPRI,'(//A/A,I5/A)') 644 & 'ERROR, analytical calculation of molecular Hessian is '// 645 & 'not available for this type of wave function.', 646 & ' Analytical derivative order available:',NAORDR, 647 & ' (You might want to choose numerical Hessian if your '// 648 & 'molecule is not too big.)' 649 CALL QUIT( 650 & 'ERROR, analytical calculation of molecular Hessian is '// 651 & 'not available for this type of wave function') 652 END IF 653 IF (MOLGRD .AND. NAORDR .LT. 1) THEN 654 WRITE(LUPRI,'(//A/A,I5)') 655 & 'ERROR, analytical calculation of molecular gradient is '// 656 & 'not available for this type of wave function.', 657 & ' Analytical derivative order available:',NAORDR 658 CALL QUIT( 659 & 'ERROR, analytical calculation of molecular gradient is '// 660 & 'not available for this type of wave function') 661 END IF 662 663C Now execute abacus 664 665 CALL ABACTL(WORK(1),LMWORK) 666 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ABACUS') 667C 668C Several files left open, have to be closed manually. 669C 670 NEWPRP = .TRUE. 671 NEWMAT = .TRUE. 672 IF (LUINTM .GT. 0) CALL GPCLOSE(LUINTM,'KEEP') 673C 674 CALL TITLER('End of Static Property Section (ABACUS) -', 675 & ' ',200) 676 CALL FLSHFO(LUPRI) 677 CALL QEXIT('ABACUS') 678 RETURN 679 END 680C /* Deck dalton_exedrv */ 681 SUBROUTINE DALTON_EXEDRV(LMWORK, WRKDLM) 682C 683C Top driver for execution of Dalton. 684C Called by PROGRAM DALTON after reading environment 685C variables. 686C 687#ifdef HAS_PCMSOLVER 688 use pcm_config, only: pcm_configuration, pcm_cfg, pcm_finalize 689#endif 690 use pelib_interface, only: use_pelib, pelib_ifc_finalize, 691 & pelib_ifc_deactivate, pelib_ifc_do_mep, 692 & pelib_ifc_do_mep_noqm, 693 & pelib_ifc_mep_noqm, pelib_ifc_do_lf 694 695#include "implicit.h" 696#include "dummy.h" 697#include "priunit.h" 698#include "mxcent.h" 699#include "maxorb.h" 700 PARAMETER (D0 = 0.0D0) 701C 702#include "gnrinf.h" 703#include "nuclei.h" 704#include "infinp.h" 705#include "inftap.h" 706#if defined (VAR_MPI) 707 INCLUDE 'mpif.h' 708#endif 709#include "infpar.h" 710#include "inftra.h" 711#include "abainf.h" 712#include "molinp.h" 713#include "cbiwlk.h" 714#include "cbirea.h" 715#include "cbisol.h" 716#include "pcmlog.h" 717Cef begin 718#include "incore.h" 719Cef end 720C 721Cholesky 722#include "ccdeco.h" 723#include "center.h" 724Cholesky 725C 726! 727! Programmers please note: 728! WORK(0) and WORK(LMOWRK+1) are allocated and contain the value 729! WRKDLM (set in main program). This is a simple (not fool-proof) test 730! that Dalton has not been out of bounds. Thus all subroutines called 731! here needs to be called with WORK(1) (and not just WORK), __except__ 732! the subroutines EXE* (here: work(0)). In these driver routines for some of the 733! major modules in Dalton, the WRKDLM is checked with BNDCHK (which is 734! in this file). /HJAaJ, Jan. 2011 735! 736 ! DIMENSION WORK(0:LMWORK) 737 738 DIMENSION CORBKP(3*MXCENT) 739 LOGICAL NEWWLK, EX, lucita_done 740 CHARACTER FILENM*12, WORD*7 741 742#include "trkoor.h" 743#include "taymol.h" 744 REAL*8, allocatable :: GRDMOL(:), HESMOL(:,:) 745 746 real(8), allocatable :: work(:) 747 748 CALL QENTER('DALTON') 749 CALL GETTIM(CSTR,WSTR) 750C 751! allocate memory using f90 utilities 752 allocate(work(0:lmwork+1),stat=i) 753! Set memory traps 754 work(0) = wrkdlm 755 work(1+lmwork) = wrkdlm 756 757 if(i /= 0)then 758 write (luerr,*) mynum,': ALLOCATE command failed to allocate'// 759 & ' the requested memory. Error code:',I 760 write (luerr,*) 'Reduce the memory demands and be welcome back' 761 call quit('Failed to allocate memory') 762 end if 763 764 IPRINT = IPRUSR 765 IPRUSR_orig = IPRUSR 766C 767 WLKREJ = .FALSE. 768 IF (OPTNEW) THEN 769 IF (QM3) 770 & CALL QUIT('Geometry optimization not implemented for .QM3') 771 RNHERM = .TRUE. 772 RNABAC = .TRUE. 773 DOWALK = .FALSE. 774 CALL OPTMIN(WORK(1),LMWORK,WRKDLM) 775 ELSE IF (OPTWLK) THEN 776 IF (QM3) 777 & CALL QUIT('ERROR: .WALK is not implemented for .QM3') 778 RNHERM = .TRUE. 779 RNABAC = .TRUE. 780 IF (NMWALK) THEN 781 CALL TITLER('Commencing numerical differentiation '// 782 & 'based on .NMDDRV',' ',200) 783 MOLGRD = .FALSE. 784 ELSE IF (IWKTYP .EQ. 6) THEN 785 CALL TITLER('Commencing numerical differentiation '// 786 & 'based on .WALK',' ',200) 787 DOWALK = .TRUE. 788 MOLGRD = .TRUE. 789 MOLHES = .TRUE. 790 ELSE 791 CALL TITLER( 792 & 'Commencing geometry optimization based on .WALK',' ',200) 793 DOWALK = .TRUE. 794 MOLGRD = .TRUE. 795 MOLHES = .TRUE. 796 END IF 797 START = .TRUE. 798 GEOCNV = .FALSE. 799 NEWWLK = .TRUE. 800 IPREAD_orig = IPREAD 801 IPREAD_reduced = -2 802C --------------------------------------- 803C loop over geometries 804C --------------------------------------- 805 10 CONTINUE 806 IF (LUCME.GT.0) REWIND (LUCME) 807 SOLVNT = .FALSE. 808 WRITE (LUPRI,'(//5X,A,I5,A//)') 809 & '@@---------- Geometry walk iteration number:',ITERNR, 810 & ' ----------' 811 CALL EXEHER(WORK(0),LMWORK,WRKDLM) 812 IPREAD = IPREAD_reduced 813C ... reduce output from READIN in following geometries 814 WLKREJ = .FALSE. 815 NEWGEO = .TRUE. 816Cef begin 817 IF (AOSAVE) THEN 818 LINTSV = .FALSE. 819 LINTMP = .FALSE. 820 INITX = .FALSE. 821 MSAVE = .TRUE. 822 MMCORE = MMWORK 823 LMCORE = MMCORE 824 ISCORE = 1 825 JSCORE = ISCORE 826 N_SHL = 1 827 I_SHL = 1 828 INDX_SHL1 = 0 829 INDX_SHL2 = 0 830 INDX_SHL3 = 0 831 INDX_SHL4 = 0 832C CALL CLEAR_INCOREMEM() 833 END IF 834Cef end 835 CALL EXESIR(WORK(0),LMWORK,WRKDLM) 836 IF (DOCCSD) CALL EXECC(WORK(0),LMWORK,WRKDLM) 837C 838C If this is a solvent run, update information for abacus 839C Add cavity center if solvent model (921014-kvm/hjaaj) 840C ===================================================== 841C 842 SOLVNT = FLAG(16) 843 IF (SOLVNT) THEN 844 NEWGEO = .TRUE. 845 NUCIND = NUCIND + 1 846 NCTOT = NCTOT + 1 847 NUCDEP = NUCDEP + 1 848 NATOMS = NATOMS + 1 849 NCNTCV = NCTOT 850 NCLINE(NUCIND) = 0 851 NAMN(NUCIND) = 'cav ' 852 NAMEX(3*NUCIND-2) = 'cav x' 853 NAMEX(3*NUCIND-1) = 'cav y' 854 NAMEX(3*NUCIND) = 'cav z' 855 NAMDEP(NUCDEP) = 'cavity' 856 NAMDPX(3*NUCDEP-2) = 'cavity x' 857 NAMDPX(3*NUCDEP-1) = 'cavity y' 858 NAMDPX(3*NUCDEP ) = 'cavity z' 859 IF (NUCDEP .GT. MXCENT) THEN 860 WRITE (LUPRI,'(//2A,/A,I5)') 861 & ' Too many atomic centers: MXCENT exceeded for', 862 & ' solvent cavity,',' Current limit:',MXCENT 863 CALL QUIT('*** ERROR *** MXCENT exceeded') 864 END IF 865 CORD(1,NUCIND) = D0 866 CORD(2,NUCIND) = D0 867 CORD(3,NUCIND) = D0 868 ISTBNU(NUCIND) = 7 869 CHARGE(NUCIND) = D0 870 CALL NUCPRO(WORK(1),LMWORK) 871CLF 872 ELSE IF(PCM) THEN 873 NPCMIN =.TRUE. 874 END IF 875C 876 IF (WLKREJ) THEN 877 WRITE (LUPRI,'(/A)') 878 & '@ Geometry step was rejected, must backstep ...' 879C 880C We have to backstep, get old geometry and call for an update 881C 882 REJECT = .TRUE. 883 CALL DCOPY(3*MXCENT,CORBKP,1,CORD,1) 884 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ',' ',IDUMMY, 885 & .FALSE.) 886 CALL WLKDRV(DUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 887 & WORK(1),LMWORK) 888 CALL GPCLOSE(LUSIFC,'KEEP') 889 REJECT = .FALSE. 890 RDINPC = .FALSE. 891 RDMLIN = .TRUE. 892C 893C Punch MOLECULE input with updated coordinates to XXXX_mol.inp 894C 895 CALL PNCMOL(ITERNR,IPRINT) 896 IF (LUSIFC .GT. 0) CALL GPCLOSE(LUSIFC,'KEEP') 897 GOTO 10 898 ELSE 899C 900C We take a backup of the present geometry in case of backsteps 901C 902 CALL DCOPY(3*MXCENT,CORD,1,CORBKP,1) 903 END IF 904C 905C We have to check if properties are to be calculated in first run 906C 907C 908C *** Walk to make numerical derivatives. *** 909C 910 IF (NMWALK) THEN 911 CALL TITLER( 912 & 'Geometry steps to make numerical derivatives', 913 & ' ',200) 914 CALL NMDINI(IPRINT) 915 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED', 916 & IDUMMY,.FALSE.) 917 REWIND (LUCMD,IOSTAT=IOS) 918 1100 READ (LUCMD,'(A7)',IOSTAT=IOS) WORD 919 IF (IOS .NE. 0) THEN 920 WRITE(LUPRI,'(//A)') 'ERROR: **NMDDRV input'// 921 & ' required when .NMDDRV requested' 922 CALL QUIT('ERROR: **NMDDRV input'// 923 & ' required when .NMDDRV requested') 924 END IF 925 CALL UPCASE(WORD) 926 IF (WORD .NE. '**NMDDR') GOTO 1100 927 CALL NMDINP(WORD,IPRINT) 928 CALL GPCLOSE(LUCMD,'KEEP') 929 CALL NUMDRV(WORK(1),LMWORK,IPRINT,WRKDLM) 930 GEOCNV = .TRUE. 931c#if defined (VAR_MPI) 932c IF (NODTOT .GE. 1 .AND. MYNUM .EQ. 0) THEN 933c IJOB = 0 934c CALL MPI_BCAST(IJOB,1,my_MPI_INTEGER,MASTER, 935c & MPI_COMM_WORLD,IERR) 936c GOTO 87 937c END IF 938c#endif 939 ELSE 940 WRINDX = .TRUE. 941 IF (NEWWLK) THEN 942 CALL WLKINI 943 IF (WALKIN) THEN 944C 945C We made some changes to walk default. Need to process these 946C 947 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ', 948 & 'FORMATTED',IDUMMY,.FALSE.) 949 REWIND (LUCMD,IOSTAT=IOS) 950 901 READ (LUCMD,'(A7)') WORD 951 CALL UPCASE(WORD) 952 IF (WORD .NE. '*WALK ') GOTO 901 953 CALL WLKINP(WORD) 954 CALL GPCLOSE(LUCMD,'KEEP') 955 END IF 956 END IF 957 IF (ITERNR .EQ. 0) THEN 958 CALL ABAINP('**START',WORK(1),LMWORK) 959 ELSE 960 CALL ABAINP('**EACH ',WORK(1),LMWORK) 961 END IF 962 IF (NEWWLK) THEN 963 NEWWLK = .FALSE. 964 CALL WLKINI 965 IF (WALKIN) THEN 966C 967C We made some changes to walk default. Need to process these 968C 969 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ', 970 & 'FORMATTED',IDUMMY,.FALSE.) 971 REWIND (LUCMD,IOSTAT=IOS) 972 900 READ (LUCMD,'(A7)') WORD 973 CALL UPCASE(WORD) 974 IF (WORD .NE. '*WALK ') GOTO 900 975 CALL WLKINP(WORD) 976 CALL GPCLOSE(LUCMD,'KEEP') 977 END IF 978 END IF 979C We may have more DOSYM under **START than under **PROPE 980 DO I = 0, 7 981 DOREPW(I) = DOSYM(I + 1) .AND. DOREPW(I) 982 END DO 983 IF (REUSED .AND. ITERNR .GT. 0) THEN 984 allocate ( GRDMOL(NCOOR), HESMOL(NCOOR,NCOOR) ) 985 CALL ABAREAD_TAYMOL(ERGMOL,GRDMOL,HESMOL,NCOOR) 986 REWIND (LUSIFC) 987 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 988 READ (LUSIFC) POTNUC,EMY,EACTIV,ERGMOL 989 CALL ABAWRIT_TAYMOL(ERGMOL,GRDMOL,HESMOL,NCOOR) 990 deallocate ( GRDMOL, HESMOL ) 991 CALL WLKDRV(DUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 992 & WORK(1),LMWORK) 993 ELSE 994 CALL EXEABA(WORK(0),LMWORK,WRKDLM) 995 END IF 996 IF (.NOT. USRIPR .AND. ((IWKTYP .EQ. 3) 997 & .OR. (IWKTYP .EQ. 6))) THEN 998 IF (MOD(ITERNR+2,10) .EQ. 0) THEN 999 IPRUSR = 0 1000 ELSE 1001 IPRUSR = -2 1002 END IF 1003 END IF 1004 WRITE (LUPRI,'(/A,I4,A,F14.6,A,F10.6)') 1005 & '@ Geometry walk step no.',ITERNR, 1006 & ': Molecular energy =', ERGMOL, 1007 & ' Molecular gradient =', GRADML 1008 END IF 1009C 1010C Punch MOLECULE input with updated coordinates to XXXX_mol.inp 1011C 1012 CALL PNCMOL(ITERNR,IPRINT) 1013C 1014 START = .FALSE. 1015 ITERNR = ITERNR + 1 1016 RDMLIN = .TRUE. 1017 IF ((.NOT. GEOCNV).AND.(ITERNR .LE. ITERMX)) THEN 1018 RDINPC = .FALSE. 1019 CALL GPINQ('RSPVEC','EXIST',EX) 1020 IF (EX) THEN 1021 LU = -9982 1022 CALL GPOPEN(LU,'RSPVEC','OLD',' ','UNFORMATTED',IDUMMY, 1023 & .FALSE.) 1024 CALL GPCLOSE(LU,'DELETE') 1025 END IF 1026 GOTO 10 1027 END IF 1028C 1029C --------------------------------------------------------------- 1030C 1031 IF (GEOCNV) THEN 1032 IPRUSR = IPRUSR_orig 1033 IF (IWKTYP .EQ. 6) THEN 1034 CALL TITLER('Nummerical differentiation finished',' ',200) 1035 ELSE 1036 CALL TITLER('Geometry walk has converged',' ',200) 1037 END IF 1038C 1039C We need to read in starting geometry before doing the vibrational 1040C analysis in a ROA calculation, K.Ruud and G.Hangartner, Oct.-96 1041C 1042 WRINDX = .TRUE. 1043 IF (IWKTYP .EQ. 6 .AND. NUMHES) THEN 1044 NEWGEO = .TRUE. 1045 CALL ABAINP('**PROPE',WORK(1),LMWORK) 1046 CALL VIBCTL(WORK(1),LMWORK) 1047 ELSE IF (IWKTYP .EQ. 6) THEN 1048 RDINPC = .FALSE. 1049 SOLVNT = .FALSE. 1050 IPREAD = IPREAD_orig 1051 CALL EXEHER(WORK(0),LMWORK,WRKDLM) 1052 NEWGEO = .TRUE. 1053 CALL EXESIR(WORK(0),LMWORK,WRKDLM) 1054 SOLVNT = FLAG(16) 1055 IF (SOLVNT) THEN 1056 NEWGEO = .TRUE. 1057 NUCIND = NUCIND + 1 1058 NUCDEP = NUCDEP + 1 1059 NATOMS = NATOMS + 1 1060 NCNTCV = NUCIND 1061 NCLINE(NUCIND) = 0 1062 NAMN(NUCIND) = 'cav ' 1063 NAMEX(3*NUCIND-2) = 'cav x' 1064 NAMEX(3*NUCIND-1) = 'cav y' 1065 NAMEX(3*NUCIND) = 'cav z' 1066 NAMDEP(NUCDEP) = 'cavity' 1067 NAMDPX(3*NUCDEP-2) = 'cavity x' 1068 NAMDPX(3*NUCDEP-1) = 'cavity y' 1069 NAMDPX(3*NUCDEP ) = 'cavity z' 1070 IF (NUCDEP .GT. MXCENT) THEN 1071 WRITE (LUPRI,'(//2A,/A,I5)') 1072 & ' Too many atomic centers: MXCENT exceeded for', 1073 & ' solvent cavity,',' Current limit:',MXCENT 1074 CALL QUIT('*** ERROR *** MXCENT exceeded') 1075 END IF 1076 CORD(1,NUCIND) = D0 1077 CORD(2,NUCIND) = D0 1078 CORD(3,NUCIND) = D0 1079 ISTBNU(NUCIND) = 7 1080 CHARGE(NUCIND) = D0 1081 CALL NUCPRO(WORK(1),LMWORK) 1082 END IF 1083 END IF 1084 NEWGEO = .TRUE. ! EXESIR resets NEWGEO to .false. 1085 ! but we need to copy abacus info to slaves 1086 CALL ABAINP('**PROPE',WORK(1),LMWORK) 1087 CALL EXEABA(WORK(0),LMWORK,WRKDLM) 1088C 1089C We also check to see if we have requested a RESPONSE calculations 1090C 1091 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED', 1092 & IDUMMY,.FALSE.) 1093 REWIND (LUCMD) 1094C 1095C Let us first see if we can find any **RESPONSE input 1096C 1097 133 CONTINUE 1098 READ (LUCMD,'(A7)',END=134,ERR=134) WORD 1099 CALL UPCASE(WORD) 1100 IF (WORD .EQ. '**RESPO') THEN 1101 CALL TITLER( 1102 & 'Starting in Dynamic Property Section (RESPONS)',' ',200) 1103 REWIND (LUCMD) 1104 RNRESP = .TRUE. 1105 CALL GPCLOSE(LUCMD,'KEEP') 1106 CALL RSPDRV(WORK(1),LMWORK) 1107 CALL TITLER( 1108 & 'End of Dynamic Property Section (RESPONS)',' ',200) 1109 ELSE 1110 GOTO 133 1111 END IF 1112 134 CONTINUE 1113 IF (LUCMD .GT. 0) CALL GPCLOSE(LUCMD,'KEEP') 1114 ELSE IF (ITERNR .GT. ITERMX) THEN 1115 CALL TITLER('WARNING - Geometry not converged',' ',200) 1116 WRITE (LUPRI,'(/A//)') '@ *** WARNING: '// 1117 & 'Maximum number of geometry iterations exceeded.' 1118 END IF 1119 ELSE 1120C ... not OPTWLK and not OPTNEW, i.e. not geometry optimization 1121 NEWGEO = .TRUE. 1122Chj IF (RNHERM .OR. TESTIN) CALL EXEHER(WORK,LMWORK,WRKDLM) 1123 CALL EXEHER(WORK(0),LMWORK,WRKDLM) 1124Chj always call EXEHER to read MOLECULE.INP (not read in abactl any more) 1125Chj instead "IF (RNHERM) CALL HERCTL" inside EXEHER /April 2009 hjaaj 1126C 1127C Request effective dipole integrals using PE library 1128 IF (USE_PELIB()) THEN 1129 IF (PELIB_IFC_DO_LF()) THEN 1130 WRITE(LUPRI,'(/A)') 'PE-QM calculation includes '// 1131 & 'effective external field effects.' 1132 END IF 1133 END IF 1134Classical MEP calculations using PE library 1135 IF (PELIB_IFC_DO_MEP() .AND. PELIB_IFC_DO_MEP_NOQM()) THEN 1136 CALL PELIB_IFC_MEP_NOQM() 1137 END IF 1138 1139Cholesky 1140 IF (CHOINT) CALL EXECHO(WORK(0),LMWORK,WRKDLM) 1141Cholesky 1142 IF (RNSIRI) CALL EXESIR(WORK(0),LMWORK,WRKDLM) 1143Cact 1144 IF (ACTSEL .AND. DOCCSD) CALL EXEACT(WORK(0),LMWORK,WRKDLM) 1145Cact 1146 IF (DOCCSD) CALL EXECC (WORK(0),LMWORK,WRKDLM) 1147C 1148C If this is a solvent run, update information for abacus 1149C Add cavity center if solvent model (921014-kvm/hjaaj) 1150C ===================================================== 1151C 1152 SOLVNT = FLAG(16) 1153 IF (SOLVNT) THEN 1154 NEWGEO = .TRUE. 1155 NUCIND = NUCIND + 1 1156 NUCDEP = NUCDEP + 1 1157 NATOMS = NATOMS + 1 1158 NCNTCV = NUCIND 1159 NCLINE(NUCIND) = 0 1160 NAMN(NUCIND) = 'cav ' 1161 NAMEX(3*NUCIND-2) = 'cav x' 1162 NAMEX(3*NUCIND-1) = 'cav y' 1163 NAMEX(3*NUCIND) = 'cav z' 1164 NAMDEP(NUCDEP) = 'cavity' 1165 NAMDPX(3*NUCDEP-2) = 'cavity x' 1166 NAMDPX(3*NUCDEP-1) = 'cavity y' 1167 NAMDPX(3*NUCDEP ) = 'cavity z' 1168 IF (NUCDEP .GT. MXCENT) THEN 1169 WRITE (LUPRI,'(//2A,/A,I5)') 1170 & ' Too many atomic centers: MXCENT exceeded for', 1171 & ' solvent cavity,',' Current limit:',MXCENT 1172 CALL QUIT('*** ERROR *** MXCENT exceeded') 1173 END IF 1174 CORD(1,NUCIND) = D0 1175 CORD(2,NUCIND) = D0 1176 CORD(3,NUCIND) = D0 1177 ISTBNU(NUCIND) = 7 1178 CHARGE(NUCIND) = D0 1179 CALL NUCPRO(WORK(1),LMWORK) 1180 END IF 1181 WRINDX = .TRUE. 1182 IF (RNRESP) THEN 1183 CALL TITLER( 1184 & 'Starting in Dynamic Property Section (RESPONS)',' ',200) 1185 CALL RSPDRV(WORK(1),LMWORK) 1186 CALL TITLER( 1187 & 'End of Dynamic Property Section (RESPONS)',' ',200) 1188 END IF 1189 IF (RNABAC) THEN 1190 CALL ABAINP('**PROPE',WORK(1),LMWORK) 1191 CALL WLKINI 1192 IF (WALKIN) THEN 1193C 1194C User has made some changes to walk default. Need to process these 1195C 1196 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED', 1197 & IDUMMY,.FALSE.) 1198 REWIND (LUCMD,IOSTAT=IOS) 1199 940 READ (LUCMD,'(A7)') WORD 1200 CALL UPCASE(WORD) 1201 IF (WORD .NE. '*WALK ') GOTO 940 1202 CALL WLKINP(WORD) 1203 CALL GPCLOSE(LUCMD,'KEEP') 1204 END IF 1205 CALL EXEABA(WORK(0),LMWORK,WRKDLM) 1206 END IF 1207 END IF 1208 87 CONTINUE 1209 1210C 1211C We finally punch information for Gamess-US graphic packages 1212C 1213 IF (SOLVNT) NUCIND = NUCIND - 1 1214 KKIND = 1 1215 KKOL = KKIND + MXCENT 1216 KLAST = KKOL + MXCENT 1217 LLEFT = LMWORK - KLAST 1218C 1219 CALL MOLPLT(WORK(KKIND),WORK(KKOL),TITMOL,WORK(KLAST),LLEFT) 1220 CALL PLTORB(WORK(KKIND),LMWORK) 1221 1222C 1223! release memory used in PElib and deactivate it 1224 IF (USE_PELIB()) THEN 1225 CALL PELIB_IFC_DEACTIVATE() 1226 CALL PELIB_IFC_FINALIZE() 1227 END IF 1228 1229 CALL GETTIM(CEND,WEND) 1230 CTOT = CEND - CSTR 1231 WTOT = WEND - WSTR 1232C 1233! Finalize pcm_scf if an EXTPCM calculation (PCM using PCMSolver 1234! library) was requested. 1235#ifdef HAS_PCMSOLVER 1236 if (pcm_cfg%do_pcm) then 1237 call pcm_finalize 1238 write(lupri, '(//A/)') 'PCMSolver interface correctly finalized' 1239 end if 1240#endif 1241 1242 CALL TIMTXT(' Total CPU time used in DALTON:',CTOT,LUPRI) 1243 CALL TIMTXT(' Total wall time used in DALTON:',WTOT,LUPRI) 1244C 1245C Stamp date and time and hostname to output 1246C 1247 CALL TSTAMP(' ',LUPRI) 1248 1249! release dynamically allocated memory 1250 deallocate(work) 1251C 1252 CALL QEXIT('DALTON') 1253 RETURN 1254 END 1255C /* Deck gnrlin */ 1256 SUBROUTINE DALTON_GNRLINP 1257C 1258C GENERAL input for Dalton 1259C 1260#ifdef HAS_PCMSOLVER 1261 use pcm_config 1262 use pcm_write 1263#endif 1264 use qmcmm_input, only: qmnpinp 1265 use pelib_interface, only: use_pelib, pelib_ifc_activate, 1266 & pelib_ifc_input_reader 1267 use qfitlib_interface, only : qfitlib_ifc_input_reader 1268 use fde_mod, only : fde_input_init, fde_dalton_input 1269 1270#include "implicit.h" 1271#include "dummy.h" 1272#include "mxcent.h" 1273#include "maxorb.h" 1274#include "priunit.h" 1275#include "gnrinf.h" 1276C 1277#include "molde.h" 1278#include "abainf.h" 1279#include "exeinf.h" 1280#include "inftra.h" 1281#include "huckel.h" 1282#include "orgcom.h" 1283#include "clsfmm.h" 1284#include "cbieri.h" 1285#include "veclen.h" 1286#include "numder.h" 1287#include "pcmdef.h" 1288#include "pcm.h" 1289#include "pcmlog.h" 1290CSONIA 1291#include "grdccpt.h" 1292#include "infpar.h" 1293C 1294Cholesky 1295#include "ccdeco.h" 1296#include "center.h" 1297C 1298 PARAMETER (NDIR = 15, NTABLE = 37, D0 = 0.D0) 1299#if defined (VAR_VECTOR) 1300 PARAMETER (IVECDF = 128) 1301#else 1302 PARAMETER (IVECDF = 1) 1303#endif 1304C 1305 LOGICAL ALLOPT 1306 1307#ifdef HAS_PCMSOLVER 1308 logical extpcm_input_provided 1309#endif 1310 1311C 1312 CHARACTER WORD*7, PROMPT*1, TABDIR(NDIR)*7, TABLE(NTABLE)*7, 1313 * REWORD*12, RWORD*6, WORD1*7 1314 DATA TABDIR/'*END OF', '*PARALL', '*WALK ', '*OPTIMI', '*PELIB ', 1315 & '*MOLBAS', '*PCM ', '*QM3 ', '*QMMM ', '*PEQM ', 1316 & '*READIN', '*QMNPMM', '*QFIT ', '*EXTPCM', '*FDE '/ 1317 DATA TABLE /'.PRINT ', '.ITERAT', '.PRIERR', 1318 & '.INTEGR', '.WAVE F', '.PROPER', ! 6 1319 & '.INPTES', '.PARALL', '.DIRECT', 1320 & '.WALK ', '.MAX IT', '.RESPON', ! 12 1321 & '.RUN PR', '.RUN RE', '.RUN WA', 1322 & '.PRESOR', '.WESTA ', '.CHOLES', ! 18 1323 & '.EXTPCM', '.FCKTRA', '.2TRATY', 1324 & '.PELIB ', '.OPTIMI', '.TOTSYM', ! 24 1325 & '.VECLEN', '.RUNERI', '.NEWTRA', 1326 & '.RUN AL', '.NMDDRV', '.PARNMD', ! 30 1327 & '.THRRED', '.DOUGLA', '.PEQM ', 1328 & '.QFIT ', '.LSLIB ', '.FDE ', ! 36 1329 & '.NOATMD'/ 1330C 1331 CALL QENTER('DALTON_GNRLINP') 1332C 1333C Initialize /VECLEN/ 1334C 1335 IVECLN = IVECDF 1336C 1337C Initialize /ABAINF/ 1338C 1339 MOLGRD = .FALSE. 1340 MOLHES = .FALSE. 1341 HELFEY = .FALSE. 1342 DOWALK = .FALSE. 1343C 1344C Initialize /MOLDEN/ 1345C 1346 MOLDEN = .TRUE. 1347C 1348C Initialize /GNRINF/ 1349C 1350 PANAS = 0.0D0 1351 THR_REDFAC = -1.0D0 ! negative number signals "not set" /hjaaj 1352 ITERNR = 0 1353 ITERMX = 20 1354 IPRUSR = 0 1355 DOFDE = .FALSE. 1356 SEGBAS = .TRUE. 1357 SEGAUX = .TRUE. 1358 WALKIN = .FALSE. 1359 HRINPC = .FALSE. 1360 SRINPC = .FALSE. 1361 RDINPC = .FALSE. 1362 SIR_INPPRC = .FALSE. 1363 RDMLIN = .FALSE. 1364 TESTIN = .FALSE. 1365 OPTWLK = .FALSE. 1366 USRIPR = .FALSE. 1367 RNHERM = .FALSE. 1368 RUNERI = .FALSE. 1369 RNSIRI = .FALSE. 1370 RNABAC = .FALSE. 1371 RNRESP = .FALSE. 1372 DOCCSD = .FALSE. 1373 WESTA = .FALSE. 1374 QMNPMM = .FALSE. 1375 NOATMD = .FALSE. 1376 1377C FCKTRA_TYPE and NEWTRA are in /INFTRA/ for practical reasons 1378C NEWTRA and NEWTRA_USEDRC are in /INFTRA/ for practical reasons 1379 FCKTRA_TYPE = -1 1380 NEWTRA = .FALSE. ! if true, use "new" transformation for Mulliken integrals 1381 NEWTRA_USEDRC = .FALSE. ! if NEWTRA and NEWTRA_USEDRC both true, also use "new" transformation for Dirac integrals 1382 1383 OPTNEW = .FALSE. 1384 NEWSYM = .FALSE. 1385 NEWBAS = .TRUE. 1386 NEWPRP = .TRUE. 1387 TOTSYM = .FALSE. 1388 RELCAL = .FALSE. 1389 NMWALK = .FALSE. 1390 DKTRAN = .FALSE. 1391 DKHINT = .FALSE. 1392 GEOALL = .FALSE. 1393 NPCMIN = .TRUE. 1394 QFIT = .FALSE. 1395 USE_LSLIB = .FALSE. 1396 IF (NODTOT .GE. 1) THEN 1397 PARCAL = .TRUE. 1398 DIRCAL = .TRUE. 1399 PARHER = .TRUE. ! in infpar.h - calculate 2-el. integrals in parallel 1400 ELSE 1401 PARCAL = .FALSE. 1402 DIRCAL = .FALSE. 1403 PARHER = .FALSE. 1404 END IF 1405C 1406C ............... 1407C Flags for calculation of long and short range integrals 1408C (mostly programmed by Jesper Kielberg Pedersen) /hjaaj 1409C For ERFEXP false: 1410C g(R12) = g_erf(R12) = erf(CHIVAL*R12) / R12 1411C where CHIVAL .ge. 0, and CHIVAL.eq.0 is normal 1/R12 1412C For ERFEXP true: 1413C g(R12) = g_erfexp(R12) 1414C = g_erf(R12) - 2*CHIVAL/SQRT(PI) * EXP(-CHIVAL**2*R12**2/3.0D0) 1415 ERFEXP(0:2) = .FALSE. 1416C 1417C IF (CHIVAL .EQ. -1.0D0) THEN calculate 1/R12 1418C ELSE calculate g(R12) listed above 1419C 1420 CHIVAL = -1.0D0 1421 CHI1ST = .FALSE. 1422C 1423C SRINTS: calculate 1/R12 - g(R12) 1424C 1425 DOSRIN = .FALSE. 1426 SRINTS = .FALSE. 1427C ............... 1428C 1429C Initialize /EXEINF/ 1430C 1431 FTWOXP = .TRUE. 1432 FTRCTL = .FALSE. 1433C ... FTRCTL true forces AO sort and integral transformation, 1434C but we may restart with old integrals. FTRCTL is set 1435C true by EXEHER signalling that new AO integrals have 1436C been calculated. 1437 NEWCMO = .TRUE. 1438 ! force integral transformation in abaset.F / abarspn.F; 1439 ! NEWCMO is then set to .false. after integral transformation. 1440 ! Will be reset to .true. each time ABACTL is called (typically: new geometry) 1441 ITRLVL_LAST = -999 ! Mulliken (Coulomb) MO integrals not available 1442 LVLDRC_LAST = -999 ! Dirac (exchange) MO integrals not available 1443C 1444C Initialize /PCMLOG/ and /EXTPCM/ 1445C 1446#ifdef HAS_PCMSOLVER 1447 extpcm_input_provided = .false. 1448#endif 1449 PCM = .FALSE. 1450 OUTFLD = .FALSE. 1451 NEWMAT = .TRUE. 1452 LOCFLD = .FALSE. 1453 NONEQ = .FALSE. 1454 NEQRSP = .FALSE. 1455Clf NPCMIN = DEFINED ABOVE..... 1456 OLDCEN = .FALSE. 1457C 1458C Initialize /HUCKEL/ 1459C 1460 DOHUCKEL = .TRUE. 1461 EWMO = .TRUE. 1462! HUCCNT = 1.75D0 1463 HUCCNT = 0.75D0 ! hjaaj: because S is ignored ... 1464 CALL IZERO(NHUCAO,8) 1465 CALL IZERO(IHUCPT,MXSHEL) 1466C 1467Cholesky 1468 CHOINT = .FALSE. 1469Cholesky 1470C 1471C Initialize /GRDCCPT/ 1472C 1473 IGRDCCPT = 1 1474 LGRDCCPT = .FALSE. 1475 1476C 1477C 1478 CALL TITLER('Output from DALTON general input processing','*',111) 1479C 1480C **** Find General input ***** 1481C 1482 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED',IDUMMY, 1483 & .FALSE.) 1484 REWIND (LUCMD,IOSTAT=IOS) 1485C ... IOSTAT to avoid program abort on some systems 1486C if reading input from a terminal 1487 J = 0 1488 900 READ (LUCMD,'(A7)',END=910,ERR=920) WORD 1489 J = J + 1 1490 CALL UPCASE(WORD) 1491 IF ((WORD .EQ. '**GENER') .OR. (WORD .EQ. '*GENERA') .OR. 1492 & (WORD .EQ. '**DALTO') .OR. (WORD .EQ. '*DALTON')) THEN 1493 GO TO 930 1494 ELSE 1495 GO TO 900 1496 END IF 1497 910 CONTINUE 1498 WRITE (LUPRI,'(A,I10,A)') 1499 & 'End of file on DALTON.INP, no **DALTON input found in', 1500 & J,' lines.' 1501 CALL QUIT('End of file on DALTON.INP, no **DALTON input found') 1502 920 CONTINUE 1503 WRITE (LUPRI,'(A,I10,A)') 1504 & 'Error reading line',J+1,' of DALTON.INP'// 1505 & ' before finding any **DALTON input.' 1506 CALL QUIT('Error reading DALTON.INP, no **DALTON input found') 1507 930 CONTINUE 1508 WORD1 = WORD 1509C 1510C ***** Process input for COMMON /GENINF/ ***** 1511C 1512C Table(01:35) : 1513C 01-05: '.PRINT ', '.ITERAT', '.PRIERR', '.INTEGR', '.WAVE F', 1514C 06-10: '.PROPER', '.INPTES', '.PARALL', '.DIRECT', '.WALK ', 1515C 11-15: '.MAX IT', '.RESPON', '.RUN PR', '.RUN RE', '.RUN WA', 1516C 16-20: '.PRESOR', '.WESTA ', '.CHOLES', '.EXTPCM', '.FCKTRA', 1517C 21-25: '.2TRATY', '.PELIB ', '.OPTIMI', '.TOTSYM', '.VECLEN', 1518C 26-30: '.RUNERI', '.NEWTRA', '.RUN AL', '.NMDDRV', '.PARNMD', 1519C 31-35: '.THRRED', '.DOUGLA', '.PEQM ', '.QFIT ', '.LSLIB ', 1520C 36-40: '.FDE ', '.NOATMD'. 1521C 1522 IPRSTAT = -10 1523 100 READ (LUCMD, '(A7)') WORD 1524 CALL UPCASE(WORD) 1525 PROMPT = WORD(1:1) 1526 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 1527 GO TO 100 1528 ELSE IF (PROMPT .EQ. '.') THEN 1529 DO 99 I = 1, NTABLE 1530 IF (TABLE(I) .EQ. WORD) THEN 1531 GO TO (101,102,103,104,105,106,107,108,109,110, 1532 & 111,112,113,114,115,116,117,118,119,120, 1533 & 121,122,123,124,125,126,127,128,129,129, 1534 & 131,132,133,134,135,136,137), I 1535 END IF 1536 99 CONTINUE 1537 IF (WORD .EQ. '.OPTION') THEN 1538 CALL PRTAB(NDIR,TABDIR, WORD1//' input keywords',LUPRI) 1539 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1540 GO TO 100 1541 END IF 1542 WRITE (LUPRI,'(/,3A,/)') 1543 & ' Keyword ',WORD,' not recognized in DALTON_GNRLINP.' 1544 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1545 CALL QUIT('Illegal keyword in DALTON_GNRLINP.') 1546 101 CONTINUE ! .PRINT 1547 READ (LUCMD,*) IPRUSR 1548 USRIPR = .TRUE. 1549 GO TO 100 1550 102 CONTINUE ! .ITERATions 1551 READ (LUCMD,*) ITERNR 1552 GO TO 100 1553 103 CONTINUE ! .PRIERR 1554 READ (LUCMD,*) IPRSTAT 1555 GO TO 100 1556 104 CONTINUE ! .INTEGRals 1557 RNHERM = .TRUE. 1558 GO TO 100 1559 105 CONTINUE ! .WAVE Function 1560 RNSIRI = .TRUE. 1561 GO TO 100 1562 106 CONTINUE ! .PROPERties 1563 RNABAC = .TRUE. 1564 GO TO 100 1565 107 CONTINUE ! .INPTESt 1566 TESTIN = .TRUE. 1567 GO TO 100 1568 108 CONTINUE ! .PARALLel calculation 1569 PARCAL = .TRUE. 1570 DIRCAL = .TRUE. 1571 PARHER = .TRUE. 1572 GO TO 100 1573 109 CONTINUE ! .DIRECT calculation (AO direct, AOTWOINT file is not generated unless needed by a module) 1574 DIRCAL = .TRUE. 1575 GO TO 100 1576 110 CONTINUE ! .WALK 1577 OPTWLK = .TRUE. 1578 DOWALK = .TRUE. 1579 GO TO 100 1580 111 CONTINUE ! .MAX ITerations ( geometry iterations ) 1581 READ (LUCMD,*) ITERMX 1582 GO TO 100 1583 112 CONTINUE ! .RESPONs 1584 RNRESP = .TRUE. 1585 GO TO 100 1586 113 CONTINUE ! .RUN PRoperties 1587 RNHERM = .TRUE. 1588 RNSIRI = .TRUE. 1589 RNABAC = .TRUE. 1590 GO TO 100 1591 114 CONTINUE ! .RUN REspons 1592 RNHERM = .TRUE. 1593 RNSIRI = .TRUE. 1594 RNRESP = .TRUE. 1595 GO TO 100 1596 115 CONTINUE ! .RUN WAve function 1597 RNHERM = .TRUE. 1598 RNSIRI = .TRUE. 1599 GO TO 100 1600 116 CONTINUE ! .PRESORt 1601 WRITE(LUPRI,'(//A/)') 1602 & 'INFO: .PRESORT deprecated, use .NEWTRA' 1603 NEWTRA = .TRUE. 1604 GO TO 100 1605 117 CONTINUE ! .WESTA 1606 WESTA = .TRUE. 1607 GO TO 100 1608 118 CONTINUE ! .CHOLESky 1609 RNHERM = .TRUE. 1610 CHOINT = .TRUE. 1611 DIRCAL = .TRUE. 1612 GO TO 100 1613 119 CONTINUE ! .EXTPCM 1614#ifdef HAS_PCMSOLVER 1615 extpcm_input_provided = .true. 1616#else 1617 WRITE(LUPRI,'(//A/)') 1618 $ 'ERROR: EXTPCM requested but PCM Module not compiled' 1619 CALL QUIT('Recompile including the PCM Module') 1620#endif 1621 goto 100 1622 120 CONTINUE ! .FCKTRA 1623 IF (FCKTRA_TYPE .LT. 0) FCKTRA_TYPE = 1 1624 NEWTRA = .TRUE. 1625 GO TO 100 1626 121 CONTINUE !.2TRATYpe 1627 ! "Secret option" for debugging/enhancing 2-electron integral transformation 1628 ! negative: old integral transformation 1629 ! zero: new integral transformation 1630 ! 1: standard .FCKTRA 1631 ! 2: .FCKTRA not parallel, even when PARHER 1632 ! 3: .FCKTRA not parallel with ONLY_J on AO2_JINT file 1633 ! 4: .FCKTRA not parallel, do not use ONLY_J on AO2_JINT file 1634 READ(LUCMD, *) FCKTRA_TYPE 1635 NEWTRA = FCKTRA_TYPE .GE. 0 1636 GO TO 100 1637 122 CONTINUE ! .PELIB 1638 CALL PELIB_IFC_ACTIVATE() 1639 GO TO 100 1640 123 CONTINUE ! .OPTIMIzation of geometry 1641 OPTNEW = .TRUE. 1642 GO TO 100 1643 124 CONTINUE ! .TOTSYM 1644 TOTSYM = .TRUE. 1645 GO TO 100 1646 125 CONTINUE ! .VECLEN 1647 READ (LUCMD,*) IVECLN 1648 GO TO 100 1649 126 CONTINUE ! .RUNERI 1650 RUNERI = .TRUE. 1651 GO TO 100 1652 127 CONTINUE ! .NEWTRA 1653 NEWTRA = .TRUE. 1654 FCKTRA_TYPE = 0 1655 GO TO 100 1656 128 CONTINUE ! .RUN ALl four modules 1657 RNHERM = .TRUE. 1658 RNSIRI = .TRUE. 1659 RNABAC = .TRUE. 1660 RNRESP = .TRUE. 1661 GO TO 100 1662 129 CONTINUE ! '.NMDDRV' and '.PARNMD' 1663 NMWALK = .TRUE. 1664 NOMOVE = .TRUE. 1665 OPTWLK = .TRUE. 1666 GOTO 100 1667 131 CONTINUE ! .THRRED 1668 READ (LUCMD,*) THR_REDFAC 1669 GOTO 100 1670 132 CONTINUE ! .DOUGLAs-kroll-hess scalar second order 1671 DKTRAN = .TRUE. 1672 GOTO 100 1673 133 CONTINUE ! .PEQM 1674 write(LUPRI, *) 'WARNING: the .PEQM option is deprecated,' 1675 & // 'please use .PELIB' 1676 CALL PELIB_IFC_ACTIVATE() 1677 GOTO 100 1678 134 CONTINUE ! .QFIT 1679#if defined (BUILD_QFITLIB) 1680 QFIT = .TRUE. 1681#else 1682 WRITE(LUPRI,*) 'ERROR for .QFIT: QFITLIB not enabled.' 1683 CALL QUIT('ERROR for .QFIT: QFITLIB not enabled.') 1684#endif 1685 GOTO 100 1686 135 CONTINUE ! .LSLIB 1687#if defined (BUILD_LSLIB) 1688 USE_LSLIB = .TRUE. 1689#else 1690 WRITE(LUPRI,*) 'ERROR for .LSLIB : LSLIB not enabled.' 1691 CALL QUIT('ERROR for .LSLIB : LSLIB not enabled.') 1692#endif 1693 GOTO 100 1694 136 CONTINUE ! .FDE 1695 DOFDE = .TRUE. 1696 GOTO 100 1697 ELSE IF (PROMPT .EQ. '*') THEN 1698 GO TO 999 1699 ELSE 1700 WRITE (LUPRI,'(/,3A,/)') ' Prompter "',PROMPT,'" illegal' 1701 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 1702 CALL QUIT('Illegal prompt in DALTON_GNRLINP.') 1703 END IF 1704 137 CONTINUE ! .NOATMD 1705 NOATMD = .TRUE. 1706 GOTO 100 1707C 1708C 1709 999 CONTINUE 1710 WRITE (LUPRI,'(1X,A)') SEPARATOR 1711 WRITE (LUPRI,'(3X,A,I5)') 'Overall default print level:',IPRUSR 1712 IF (IPRSTAT .EQ. -10) IPRSTAT = IPRUSR + 1 1713 WRITE (LUPRI,'(3X,A,I5/)') 'Print level for DALTON.STAT:',IPRSTAT 1714! WRITE (LUERR,'(A,I5)') ' Print level for DALTON.STAT',IPRSTAT 1715 WRITE (LUSTAT,'(A,I5)') ' Print level for DALTON.STAT',IPRSTAT 1716C ... default print level for statistics /hjaaj apr 2000 1717 IF (TESTIN) WRITE (LUPRI,'(/A/)') ' *** Input test run only ***' 1718 IF (THR_REDFAC .GT. 0.0D0) THEN 1719 WRITE (LUPRI,'(A,1P,D10.2)') 1720 & '@ INFO: All thresholds are multiplied by',THR_REDFAC 1721 IF (THR_REDFAC .LT. 1.D-10) THEN 1722 WRITE (LUPRI,'(A)') ' FATAL INPUT ERROR: Illegal value' 1723 CALL QUIT('Illegal input value for .THRRED') 1724 END IF 1725 END IF 1726 IF (PARCAL .AND. NODTOT .EQ. 0) THEN 1727 NINFO = NINFO + 1 1728 WRITE (LUPRI,'(4X,A)') SEPARATOR,'INFO: '// 1729 & 'Request for parallel calculation (.PARALL) is ignored '// 1730 & 'because only one CPU process available.',SEPARATOR 1731 PARCAL = .FALSE. 1732 PARHER = .FALSE. 1733 END IF 1734#if defined (VAR_MPI) 1735 IF (PARCAL) WRITE (LUPRI,'(4X,A)')'Parallel calculation using MPI' 1736#else 1737 IF (PARCAL) THEN 1738 WRITE (LUPRI,'(//4X,A/4X,A)') 1739 & 'Parallel calculation requested, but this is not a', 1740 & 'parallel version. Please recompile for MPI' 1741 CALL QUIT('Parallel calculation requested w/o MPI activated') 1742 END IF 1743#endif 1744 IF (OPTWLK) THEN 1745 IF (NMWALK) THEN 1746 WRITE (LUPRI,'(A/A,I5)') 1747 & '@ Numerical derivatives will be calculated.', 1748 & ' Number of different geometries:', ITERMX 1749 ELSE 1750 WRITE (LUPRI,'(A/6X,A,I5/6X,A,I5)') 1751 & '@ Geometry optimization', 1752 & 'Starting at iteration:', ITERNR, 1753 & 'Maximum number of steps:', ITERMX 1754 END IF 1755 ELSE IF (ITERNR .GT. 0) THEN 1756 NINFO = NINFO + 1 1757 WRITE (LUPRI,'(4X,A/4X,A/6X,A/4X,A)') SEPARATOR, 1758 & 'INFO: Input specification of geometry iteration number '// 1759 & 'is only valid for ".WALK ", ".NMDDRV", or ".PARNMD"', 1760 & 'Geometry iteration number has been reset to 0',SEPARATOR 1761 ITERNR = 0 1762 END IF 1763 IF (TOTSYM .AND. OPTNEW) THEN 1764 WRITE(LUPRI,'(//A/4X,A/6X,A)') SEPARATOR, 1765 & 'ERROR: .TOTSYM is incompatible with .OPTIMIZE', 1766 & 'Use instead the .WALK module together with .TOTSYM' 1767 CALL QUIT('.TOTSYM and .OPTIMIZE is incompatible options') 1768 END IF 1769 IF (TOTSYM) WRITE (LUPRI,'(4X,A)') 1770 & 'Only totally symmetric part of molecular Hessian calculated' 1771 IF (DIRCAL) WRITE(LUPRI,'(4X,A)') 1772 & 'AO-direct calculation (in sections where implemented)' 1773 IF (DKTRAN) WRITE (LUPRI,'(4X,A)') 1774 & 'The second order Douglas-Kroll-Hess transformation is applied' 1775 IF (RNHERM) WRITE (LUPRI,'(4X,A)') 1776 & 'HERMIT 1- and 2-electron integral sections will be executed' 1777 IF (RUNERI) THEN 1778 WRITE (LUPRI,'(4X,A)') 1779 & '2-elctron integrals are calculated using ERI '// 1780 & 'instead of HERMIT (where possible)' 1781! RUNERI = RUNERI .AND. .NOT. PARCAL 1782! IF (.NOT. RUNERI) THEN 1783! NINFO = NINFO + 1 1784! WRITE (LUPRI,'(3X,A)') 1785! & ' INFO: NO! ERI is only implemented for direct and not parallel' 1786! END IF 1787 END IF 1788 IF (IVECLN .NE. 1) WRITE (LUPRI,'(4X,A,I5)') 1789 & 'Vector length used in direct Fock matrix calculations '// 1790 & '(assuming vector machine)',IVECLN 1791 IF (FCKTRA_TYPE .gt. 0) THEN 1792 WRITE(LUPRI,'(4X,A,I5)') 1793 & 'Fock-matrix based 2-el integral transformation (2016) used,' 1794 & //' type',FCKTRA_TYPE 1795 ELSE IF (NEWTRA) THEN 1796 WRITE(LUPRI,'(4X,A)') 1797 & '"New" (from 1988!) integral transformation used' 1798 ELSE 1799 WRITE(LUPRI,'(4X,A)') '"Old" integral transformation used'// 1800 & ' (limited to max 255 basis functions)' 1801 END IF 1802 IF (RNSIRI) WRITE (LUPRI,'(4X,A)') 1803 & 'Wave function sections will be executed (SIRIUS module)' 1804 IF (RNRESP) WRITE (LUPRI,'(4X,A)') 1805 & 'Dynamic molecular response properties section will be'// 1806 & ' executed (RESPONSE module)' 1807 IF (RNABAC) WRITE (LUPRI,'(4X,A)') 1808 & 'Static molecular property section will be executed'// 1809 & ' (ABACUS module)' 1810 IF (NMWALK) WRITE (LUPRI,'(4X,A)') 1811 & 'Numerical derivatives will be calculated' 1812 IF (WESTA) WRITE (LUPRI,'(4X,A)') 1813 & 'Information for WESTA will be calculated and written to files.' 1814 IF (USE_PELIB()) WRITE(LUPRI,'(4X,A)') 1815 & 'Environment effects are included (PElib)' 1816 IF(QFIT) WRITE(LUPRI,'(4X,A)') 1817 & 'Potential derived multipole moments will be calculated'// 1818 & ' (QFITLIB)' 1819C 1820Cholesky 1821 IF (CHOINT) THEN 1822 IF (PARCAL) CALL PARQUIT('CHOLESKY') 1823 WRITE (LUPRI,'(4X,A)') 1824 & 'Cholesky decomposition-based calculation will be done' 1825 ENDIF 1826Cholesky 1827C 1828 WRITE (LUPRI,'(1X,A)') SEPARATOR 1829C 1830C TABDIR : '*END OF', '*PARALL', '*WALK ', '*OPTIMI', '*PELIB ', 1831C '*MOLBAS', '*PCM ', '*QM3 ', '*QMMM ', '*PEQM ', 1832C '*READIN', '*QMNPMM', '*QFIT ', '*EXTPCM', '*FDE ' 1833C 1834C Initialize for REAINP and READIN (reading of .mol input file) 1835 IPREAD_ini = IPRUSR + 1 1836 CALL REAINI(IPREAD_ini,RELCAL,TSTINP) 1837C 1838 200 PROMPT = WORD(1:1) 1839 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') THEN 1840 READ (LUCMD,'(A7)') WORD 1841 CALL UPCASE(WORD) 1842 GO TO 200 1843 ELSE IF (PROMPT .EQ. '*') THEN 1844 DO 210 I = 1, NDIR 1845 IF (WORD .EQ. TABDIR(I)) THEN 1846 GOTO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15), I 1847 END IF 1848 210 CONTINUE 1849 IF (WORD(1:2) .EQ. '**') GO TO 1 1850 WRITE (LUPRI,'(/,3A,/)') ' Directory ',WORD,' nonexistent.' 1851 CALL PRTAB(NDIR,TABDIR,WORD1//' input keywords',LUPRI) 1852 CALL QUIT('Illegal directory in DALTON_GNRLINP.') 1853 ELSE 1854 WRITE (LUPRI,'(/,3A,/)') ' Prompter "',PROMPT,'" illegal or', 1855 * ' out of order.' 1856 CALL PRTAB(NDIR,TABDIR,WORD1//' input keywords',LUPRI) 1857 CALL QUIT('Program stopped in DALTON_GNRLINP, illegal prompt.') 1858 END IF 1859C "*PARALL" input 1860 2 CONTINUE 1861#if defined (VAR_MPI) 1862 CALL PRLINP(WORD) 1863#else 1864 WRITE (LUPRI,'(//4X,A/4X,A)') 1865 & '"*PARALL" input sprecified, but this is not a parallel', 1866 & 'version. Please recompile DALTON for MPI.' 1867 CALL QUIT('*PARALL input specified w/o MPI activated') 1868#endif 1869 GO TO 200 1870C "*WALK " input 1871 3 WALKIN = .TRUE. 1872C 1873C We have to wait with the input processing of "*WALK " 1874C until we know more about the molecule 1875C 1876 203 READ (LUCMD,'(A7)') WORD 1877 IF (WORD(1:1) .NE. '*') GOTO 203 1878 GO TO 200 1879 1880C "*OPTIMI" input 1881 4 CONTINUE 1882 CALL OPINPU(WORD) 1883 GO TO 200 1884 1885C "*PELIB " input 1886 5 CONTINUE 1887 CALL PELIB_IFC_INPUT_READER(WORD) 1888 GO TO 200 1889 1890C "*MOLBAS" input 1891 6 CONTINUE 1892 CALL REAINP(WORD,RELCAL) 1893 GO TO 200 1894 1895Clf: "*PCM " input processing 1896 7 CONTINUE 1897C INPERR is defined as a precaution (not used during PCM input processing) 1898 INPERR = 0 1899 ALLOPT = .FALSE. 1900 CALL PCMINP(WORD,INPERR,ALLOPT) 1901 GO TO 200 1902C 1903C "*QM3 " input 1904 8 CONTINUE 1905 CALL QM3INP(WORD) 1906 GO TO 200 1907C 1908C "*QMMM " input 1909 9 CONTINUE 1910 CALL QMMMIN(WORD) 1911 GO TO 200 1912 1913! "*PEQM " input 1914 10 CONTINUE 1915 write(LUPRI, *) 'WARNING: the *PEQM section is deprecated,' 1916 & // 'please use *PELIB' 1917 CALL PELIB_IFC_INPUT_READER(WORD) 1918 GO TO 200 1919C 1920C "*READIN" input ( DEPRECATED !!! Use "*MOLBAS") 1921 11 CONTINUE 1922 WRITE(LUPRI,'(//A/)') 1923 & ' INFO : "*READIN" is deprecated, use instead "*MOLBAS"' 1924 GO TO 6 1925C 1926C "*QMNPMM" input 1927 12 CONTINUE 1928 CALL QMNPINP(WORD) 1929 GO TO 200 1930 1931! "*QFIT " input 1932 13 CONTINUE 1933 CALL QFITLIB_IFC_INPUT_READER(WORD) 1934 GOTO 200 1935 1936! *EXTPCM input 1937 14 CONTINUE 1938#ifdef HAS_PCMSOLVER 1939 call move_to_next_star(word,lucmd) 1940#else 1941 WRITE(LUPRI,'(//A/)') 1942 $ 'ERROR: EXTPCM requested but PCM Module not compiled' 1943 CALL QUIT('Recompile including the PCM Module') 1944#endif 1945 GOTO 200 1946 1947! *FDE input 1948 15 CONTINUE 1949 call fde_input_init(LUPRI,'DALTON.OUT') 1950 call fde_dalton_input(word,.true.) 1951 GOTO 200 1952 1953C "*END OF" or "**something" 1954 1 CONTINUE 1955C 1956#ifdef HAS_PCMSOLVER 1957! read EXTPCM input if applicable 1958 if(extpcm_input_provided)then 1959 call read_menu_input('DALTON.INP', lucmd, 1960 & '*EXTPCM', extpcm_input_provided) 1961 pcm_cfg%do_pcm = .true. 1962 if (.not. extpcm_input_provided) then 1963 call quit('Oops, read_menu_input could not find *EXTPCM input') 1964 end if 1965 end if 1966#endif 1967 1968#if defined (VAR_MPI) 1969 CALL PRLINP(WORD) 1970#endif 1971 CALL OPINPU(WORD) 1972C 1973C In order to see if there is any change in starting orbitals (i.e. 1974C not Huckel guess, we have to search for *ORBITA 1975C 1976 REWIND (LUCMD,IOSTAT=IOS) 1977C ... IOSTAT to avoid program abort on some systems 1978C if reading input from a terminal 1979 5000 READ (LUCMD,'(A7)',END=8000) WORD 1980 CALL UPCASE(WORD) 1981 IF (WORD .EQ. '*ORBITA') THEN 1982 400 READ (LUCMD,'(A7)') WORD 1983 CALL UPCASE(WORD) 1984 IF (WORD .EQ. '.MOSTAR') THEN 1985 READ (LUCMD,'(A7)') WORD 1986 CALL UPCASE(WORD) 1987 IF (INDEX(WORD,'EWMO') .GT. 0) THEN 1988 EWMO = .TRUE. 1989 DOHUCKEL = .TRUE. 1990 ELSE IF (INDEX(WORD,'HUCKEL') .GT. 0 1991 & .OR. INDEX(WORD,'EHT') .GT. 0) THEN 1992 EWMO = .FALSE. 1993 DOHUCKEL = .TRUE. 1994 ELSE 1995 DOHUCKEL = .FALSE. 1996 END IF 1997 ELSE IF (WORD(1:1) .EQ. '*') THEN 1998 GOTO 8000 1999 ELSE 2000 GOTO 400 2001 END IF 2002 ELSE 2003 GOTO 5000 2004 END IF 2005 8000 CALL FLSHFO(LUPRI) 2006C 2007C *** Doing just a survey of the wavefunctions we use. *** 2008C 2009c INPERR = 0 2010c NUMRUN = 0 2011c IRDMO = 8 2012c IREST = 0 2013c NSYM = 1 2014c REWIND (LUCMD,IOSTAT=IOS) 2015c 8100 READ (LUCMD,'(A7)') WORD 2016c IF (.NOT.((WORD.EQ.'**WAVE ').OR.(WORD.EQ.'**WAVEF') 2017c & .OR.(WORD.EQ.'**SIRIUS'))) GOTO 8100 2018c CALL NEWINP(INPERR,NUMRUN,IRDMO,IREST,NSYM,NORB) 2019C 2020 CALL GPCLOSE(LUCMD,'KEEP') 2021 CALL QEXIT('DALTON_GNRLINP') 2022 RETURN 2023 END 2024#if defined (VAR_PARIO) 2025C /* Deck pariot */ 2026 SUBROUTINE PARIOT 2027C 2028C Master routine for checking whether we will do parallell I/O by 2029C parallelizing over nuclear geometries 2030C 2031#include "implicit.h" 2032 INCLUDE 'mpif.h' 2033#include "dummy.h" 2034#include "mxcent.h" 2035#include "maxorb.h" 2036 CHARACTER WORD*7 2037#include "priunit.h" 2038#include "inftap.h" 2039#include "molinp.h" 2040#include "infpar.h" 2041#include "pario.h" 2042C 2043 PARIO = .FALSE. 2044 CALL GPOPEN(LUCMD,'DALTON.INP','OLD',' ','FORMATTED',IDUMMY, 2045 & .FALSE.) 2046 5000 READ (LUCMD,'(A7)',END=4070) WORD 2047 CALL UPCASE(WORD) 2048 IF (WORD .EQ. '.PARNMD') THEN 2049 PARIO = .TRUE. 2050 GOTO 4070 2051 ELSE 2052 GOTO 5000 2053 END IF 2054 4070 CONTINUE 2055C 2056C We are going to do parallel I/O. Inform the slaves about this 2057C 2058 CALL MPI_BCAST(PARIO,1,my_MPI_LOGICAL,MASTER,MPI_COMM_WORLD, 2059 & IERR) 2060C 2061C Read in Dalton input. Save it temporarily MLINE 2062C 2063 IF (PARIO) THEN 2064#if defined (VAR_MPI2) 2065C 2066C The following piece of code should be replaced by an initialization 2067C of the counter for the numerical differentiation which will be addressed 2068C using RMA operations. It should probably be in a common block 2069C 2070 I0 = 1 2071 LUNMCL = -9056 2072 CALL GETENV('WRKDIR',WRKDIR) 2073 LENWRK = 0 2074 DO 43 I = 1, 60 2075 IF (WRKDIR(I:I) .EQ. ' ') GO TO 44 2076 LENWRK = LENWRK + 1 2077 43 CONTINUE 2078 44 CONTINUE 2079 IF (WRKDIR(LENWRK:LENWRK) .NE. '/') THEN 2080 LENWRK = LENWRK + 1 2081 WRKDIR(LENWRK:LENWRK) = '/' 2082 END IF 2083 CALL GPOPEN(LUNMCL,WRKDIR(1:LENWRK)//'NUMCAL','NEW',' ', 2084 & 'FORMATTED',IDUMMY,.FALSE.) 2085 WRITE (LUNMCL,'(I5)') I0 2086 CALL GPCLOSE(LUNMCL,'KEEP') 2087#endif 2088 DO ILINE = 1, KMLINE 2089 MLINE(ILINE) = ' '// 2090 & ' ' 2091 END DO 2092 REWIND (LUCMD) 2093 JLINE = 1 2094 5030 READ (LUCMD,'(A8)',END=5010) MLINE(JLINE) 2095 JLINE = JLINE + 1 2096 GOTO 5030 2097 5010 CALL MPI_BCAST(JLINE,1,my_MPI_INTEGER,MASTER, 2098 & MPI_COMM_WORLD,IERR) 2099 CALL MPI_BCAST(MLINE,80*JLINE,MPI_CHARACTER,MASTER, 2100 & MPI_COMM_WORLD,IERR) 2101 CALL GPCLOSE(LUCMD,'KEEP') 2102C 2103C Now we read and send MOLECULE.INP 2104C 2105 DO ILINE = 1, JLINE 2106 MLINE(ILINE)(1:8) = ' ' 2107 END DO 2108 CALL GPOPEN(LUMOL,'MOLECULE.INP',' ',' ','FORMATTED',IDUMMY, 2109 & .FALSE.) 2110 JLINE = 1 2111 6030 READ (LUMOL,'(A80)',END=6010) MLINE(JLINE) 2112 JLINE = JLINE + 1 2113 GOTO 6030 2114 6010 CALL MPI_BCAST(JLINE,1,my_MPI_INTEGER,MASTER, 2115 & MPI_COMM_WORLD,IERR) 2116 CALL MPI_BCAST(MLINE,80*JLINE,MPI_CHARACTER,MASTER, 2117 & MPI_COMM_WORLD,IERR) 2118 CALL GPCLOSE(LUMOL,'KEEP') 2119 ELSE 2120 REWIND (LUCMD) 2121 END IF 2122C 2123C We should have sent everything we need know. Return to the 2124C main program. 2125C 2126 RETURN 2127 END 2128C /* Deck parion */ 2129 SUBROUTINE PARION 2130C 2131C Slave routine for learning whether we are going to do parallell I/O 2132C by parallelizing over nuclear geometries 2133C 2134#include "implicit.h" 2135#include "dummy.h" 2136#include "mxcent.h" 2137#include "maxorb.h" 2138#include "priunit.h" 2139#include "inftap.h" 2140#include "molinp.h" 2141#include "infpar.h" 2142#include "pario.h" 2143 2144 INCLUDE 'mpif.h' 2145C 2146C We start by waiting for the master to tell us whether we will do 2147C parallel I/O or not 2148C 2149#ifdef VAR_INT64 2150 call quit( 2151 &'Routine PARION is not revised for int64 code and int32 mpi') 2152#endif 2153 CALL MPI_BCAST(PARIO,1,my_MPI_LOGICAL,MASTER,MPI_COMM_WORLD, 2154 & IERR) 2155 IF (PARIO) THEN 2156 CALL GETENV('WRKDIR',WRKDIR) 2157 LENWRK = 0 2158 DO 43 I = 1, 60 2159 IF (WRKDIR(I:I) .EQ. ' ') GO TO 44 2160 LENWRK = LENWRK + 1 2161 43 CONTINUE 2162 44 CONTINUE 2163 IF (WRKDIR(LENWRK:LENWRK) .NE. '/') THEN 2164 LENWRK = LENWRK + 1 2165 WRKDIR(LENWRK:LENWRK) = '/' 2166 END IF 2167 CALL GPOPEN(LUCMD,'DALTON.INP','NEW',' ','FORMATTED',IDUMMY, 2168 & .FALSE.) 2169 CALL GPOPEN(LUMOL,'MOLECULE.INP','NEW',' ','FORMATTED',IDUMMY, 2170 & .FALSE.) 2171 CALL MPI_BCAST(JLINE,1,my_MPI_INTEGER,MASTER, 2172 & MPI_COMM_WORLD,IERR) 2173 CALL MPI_BCAST(MLINE,80*JLINE,MPI_CHARACTER,MASTER, 2174 & MPI_COMM_WORLD,IERR) 2175 REWIND (LUCMD) 2176 DO ILINE = 1, JLINE 2177 WRITE (LUCMD,'(A80)') MLINE(ILINE) 2178 END DO 2179 CALL GPCLOSE(LUCMD,'KEEP') 2180C 2181C We have written DALTON.INP. Now we collect and write MOLECULE.INP 2182C 2183 CALL MPI_BCAST(JLINE,1,my_MPI_INTEGER,MASTER, 2184 & MPI_COMM_WORLD,IERR) 2185 CALL MPI_BCAST(MLINE,80*JLINE,MPI_CHARACTER,MASTER, 2186 & MPI_COMM_WORLD,IERR) 2187 REWIND (LUMOL) 2188 DO ILINE = 1, JLINE 2189 WRITE (LUMOL,'(A80)') MLINE(ILINE) 2190 END DO 2191 CALL GPCLOSE(LUMOL,'KEEP') 2192 END IF 2193 RETURN 2194 END 2195#endif /* VAR_PARIO WARNING: DO NOT activate VAR_PARIO without fixing the PARIO code */ 2196C 2197Cholesky 2198C /* Deck execho */ 2199 SUBROUTINE EXECHO(WORK, LMWORK, WRKDLM) 2200#include "implicit.h" 2201#include "dummy.h" 2202 DIMENSION WORK(0:LMWORK) 2203#include "priunit.h" 2204#include "gnrinf.h" 2205#include "dftcom.h" 2206C 2207C Run Cholesky stuff 2208C 2209 CALL QENTER('CHOLESKY') 2210C 2211 CALL TITLER('Starting in Cholesky Integral Section', 2212 & ' ',200) 2213 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini CHOLESKY') 2214C 2215 CALL CC_CHODRV(WORK(1),LMWORK) 2216 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'CHOLESKY') 2217C 2218ctst CLOSE(24) 2219 CALL TITLER('End of Cholesky Integral Section',' ',200) 2220C 2221 CALL FLSHFO(LUPRI) 2222 CALL QEXIT('CHOLESKY') 2223 RETURN 2224 END 2225C 2226C /* Deck exeact */ 2227C 2228 SUBROUTINE EXEACT(WORK, LMWORK, WRKDLM) 2229#include "implicit.h" 2230#include "dummy.h" 2231 DIMENSION WORK(0:LMWORK) 2232#include "priunit.h" 2233#include "gnrinf.h" 2234#include "mxcent.h" 2235#include "maxorb.h" 2236#include "center.h" 2237casm 2238#include "ccorb.h" 2239C 2240C 2241C Select orbitals and basis fucntions for hybrid CC models. 2242C 2243 CALL QENTER('EXEACT') 2244C 2245 CALL TITLER('Starting in Choact Section',' ',200) 2246 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'ini CHOACT') 2247C 2248 CALL CC_SELACT(WORK(1),LMWORK) 2249 CALL BNDCHK(WORK(0), LMWORK+2, WRKDLM, 'CHOACT') 2250C 2251ctst CLOSE(24) 2252 CALL TITLER('End of Choact Section',' ',200) 2253C 2254 CALL FLSHFO(LUPRI) 2255 CALL QEXIT('EXEACT') 2256 RETURN 2257 END 2258C --- end of dalgnr.F --- 2259