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! 18! 19! File: dalton/rsp/lagrang.F 20! 21! The routines in this file are generally related to 22! setup the Lagrangian for doublet reference states with 23! triplet property operators, and solving for the corresponding 24! response vector. 25! 26! This is what has been called the "restricted-unrestricted" approach 27! iin the literature. 28! 29! The driver RSPESR uses this for calculating expectation vaules of triplet 30! operators for non-singlet reference states. 31! This file also includes the input routine ESRINP for specification 32! of ESR calculations, using RSPESR here for the requested expectation values, 33! as well as routines in rspg.F and in rspzfs.F for g-tensor and ZFS, 34! depending on input. 35! 36! 37C=========================================================================== 38CRevision 1.4 2000/05/24 18:45:09 hjj 39Cnew getref calls with appropriate NDREF (fixing CSF and triplet) 40Cstop if solvent and CSF 41Cbugfix: too little was allocated for KGP2 42C=========================================================================== 43C 44C 45C /* Deck lagran */ 46 SUBROUTINE LAGRAN(WORD,FC,FV,CMO,UDV,PV,XINDX,WRK,LWRK) 47C 48C 18-AUG-1991 49C 50C Called from GETGPV when label(5:8) = 'LAGR' 51C 52#include "implicit.h" 53#include "dummy.h" 54 DIMENSION CMO(*),UDV(NASHT,*),FC(*),FV(*),PV(*),XINDX(*),WRK(*) 55 CHARACTER*8 WORD 56C-- common blocks: 57#include "maxorb.h" 58#include "priunit.h" 59#include "infdim.h" 60#include "infvar.h" 61#include "inforb.h" 62#include "infrsp.h" 63#include "wrkrsp.h" 64#include "infpri.h" 65#include "infinp.h" 66#include "dftcom.h" 67#include "pcmlog.h" 68C 69 PARAMETER (DHALF=0.5D0,D1=1.0D0,DM1=-1.0D0) 70C 71 CALL QENTER('LAGRAN') 72C 73 IF (IPRRSP .GE. 10) WRITE (LUPRI,7010) WORD 74 7010 FORMAT(//' Output from LAGRAN:'/' ===================' 75 * //' Property label = ',A/) 76 IF (SOPPA) THEN 77 WRITE (LUPRI,*) 'SOPPA not implemented in LAGRAN yet, sorry!' 78 CALL QUIT('LAGRAN-ERROR: SOPPA not implemented') 79 END IF 80clf IF ((FLAG(16).OR.PCM).AND.(NASHT.GT.1).AND.(NCREF.NE.KZCONF)) THEN 81Chj: see comments in RSPSLV in rspsol.F 82clf WRITE(LUPRI,*)'Solvent not implemented in LAGRAN with CSFs' 83clf WRITE(LUPRI,*)'Use .DETERMINANTS and try again!' 84clf CALL QUIT('Solvent not implemented in LAGRAN with CSFs') 85clf END IF 86C 87C ALLOCATE WORK SPACE 88C 89 KGP = 1 90 KCREF = KGP + KZYVAR 91 NDREF = MAX(KZCONF,NCREF) 92Chj ... RSPOLI requires determinants for ZYCVEC when triplet, 93Chj ... thus KZCONF. However, for ROHF NCREF = 1, KZCONF = 0 94Chj and we need CREF(1) for RSPDM. Thus MAX(KZCONF,NCREF) 95 KTOT = KCREF + NDREF 96C 97 CALL DZERO(WRK(KGP),KZYVAR) 98 CALL GETREF(WRK(KCREF),NDREF) 99C 100C CONSTRUCT the Lagrangian vector, the GP vector, 101C using routine for the ORBITAL PART OF LINEAR TRANSFORMED VECTORS 102C 103 IF (.NOT.RSPCI) THEN 104C 105C ALLOCATE WORK SPACE FOR RSPOLI 106C NCSIM = 1, NOSIM = 0 107C 108 KFVTD = KTOT 109 LFVTD = N2ORBX 110 IF (DOMCSRDFT) LFVTD = 2*LFVTD 111C ... Extra allocation for "FCTD" in MCSCF-SRDFT, 112C is needed for the VxcTD matrix. 113 KQATD = KFVTD + LFVTD 114 KQBTD = KQATD + NORBT * NASHT 115 KWRK1 = KQBTD + NORBT * NASHT 116 LWRK1 = LWRK - KWRK1 117 IF (LWRK1.LT.0) CALL ERRWRK('LAGRAN',KWRK1-1,LWRK) 118 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 119 CALL QUIT('LAGRANGE is not implemented yet for srDFT') 120 END IF 121 IF ((DFTADD.OR.HSROHF).AND.TRPLET) THEN 122 CALL DAXPY(NNORBT,D1,FC,1,FV,1) 123 CALL DSCAL(NNORBT,DM1,FV,1) 124 END IF 125 CALL RSPOLI(1,0,UDV,WRK(KCREF),NDREF,FC,FV,PV,DUMMY, 126 * DUMMY,DUMMY, DUMMY,DUMMY,WRK(KFVTD), 127 * WRK(KQATD),WRK(KQBTD),WRK(KGP), 128 * XINDX,CMO,WRK(KWRK1),LWRK1) 129C CALL RSPOLI(NCSIM,NOSIM,UDV,ZYCVEC,LZYCVEC,FC,FV,PVX,ZYMAT, 130C * FCONE,FVONE,QAONE,QBONE,FVTD,QATD,QBTD,EVECS, 131C * XINDX,CMO,WRK,LWRK) 132 IF ((DFTADD.OR.HSROHF).AND.TRPLET) THEN 133 CALL DAXPY(NNORBT,D1,FC,1,FV,1) 134 CALL DSCAL(NNORBT,DM1,FV,1) 135 END IF 136C 137 END IF 138 IF (IPRRSP.GT.120) THEN 139 WRITE(LUPRI,'(/A)') ' LAGRAN: LAGRANGIAN VECTOR ' 140 CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 141 END IF 142C 143C section for calculating esr-properties with solvent 144C contributions. 145C 146 IF (FLAG(16).OR.PCM) THEN 147C 148C workspace allocation 149C 150 KGP2 = KTOT 151 KTDV = KGP2 + NVARH 152Chj KTDV = KGP2 + KZVAR 153Chj ... SOLGDT uses NVARH (.gt. KZVAR !) for GP2 154 KDV = KTDV + N2ASHX 155 KWRK1 = KDV + NNASHX 156 IF (FLAG(16)) THEN 157 KWRK2 = KWRK1 + NLMSOL * 2 158 ELSE IF (PCM) THEN 159 KWRK2 = KWRK1 160 ELSE 161 KWRK2 = KWRK1 162 END IF 163 LWRK1 = LWRK - KWRK2 164C 165C initialize solvent gradient. 166C 167 CALL DZERO(WRK(KGP2),NVARH) 168C 169C 170 ISPIN1 = 1 171 ISPIN2 = 0 172 CALL RSPDM(IREFSY,IREFSY,NDREF,NDREF,WRK(KCREF),WRK(KCREF), 173 * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 174 * XINDX,WRK,KWRK2,LWRK1) 175C CALL RSPDM(ISYM,ISYM,NCDIM.NCDIM,REF,REF, 176C TUDV,DUMMY, ISPIN1,ISPIN2,TDM,NORHO2, 177C XINDX,WRK,KFREE,LFREE) 178 IF (IPRRSP.GT.110) THEN 179 WRITE(LUPRI,'(/A)')'LAGRANG: ONE ELECTRON SPIN DENSITY' 180 CALL OUTPUT(WRK(KTDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 181 END IF 182C 183C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX 184C 185 CALL DGETSP(NASHT,WRK(KTDV),WRK(KDV)) 186C 187 IF (IPRRSP.GT.50 .AND. NASHT.GT.0) THEN 188 WRITE(LUPRI,'(/A)') 189 & 'LAGRANG solvent: PACKED ONE ELECTRON SPIN DENSITY MATRIX' 190 CALL OUTPAK(WRK(KDV),NASHT,1,LUPRI) 191 END IF 192C 193C use wavefunction routine to calculate solvent gradient. 194C 195 IF (FLAG(16)) THEN 196 CALL SOLGDT(WRK(KCREF),CMO,XINDX,WRK(KDV),WRK(KGP2), 197 * ESOLT,WRK(KWRK1),WRK(KWRK2),LWRK1,NDREF,IREFSY) 198 ELSE IF (PCM) THEN 199 CALL PCMGDT(WRK(KCREF),CMO,XINDX,WRK(KDV),WRK(KGP2), 200 * ESOLT,WRK(KWRK2),LWRK1,NDREF,IREFSY) 201 ELSE 202 CALL QUIT("Lagran: Either FLAG(16) or PCM must be TRUE") 203 END IF 204 205C 206C scale by a half for consistency with response and add to GP. 207C 208 CALL DAXPY(KZWOPT,DHALF,WRK(KGP2+KZCONF),1,WRK(KGP+KZCONF),1) 209 CALL DAXPY(KZWOPT,(-DHALF),WRK(KGP2+KZCONF),1, 210 * WRK(KGP+KZVAR+KZCONF),1) 211C 212 IF (IPRRSP.GT.120) THEN 213 WRITE(LUPRI,'(/A)') ' LAGRAN: 2 * SOLVENT LAGRANGIAN VECTOR' 214 CALL OUTPUT(WRK(KGP2),1,NVARH,1,1,NVARH,1,1,LUPRI) 215 WRITE(LUPRI,'(/A)') 216 & ' LAGRAN: Electronic + solvent LAGRANGIAN VECTOR' 217 CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 218 END IF 219C 220 ENDIF 221C 222C *** END OF LAGRAN 223C 224 CALL QEXIT('LAGRAN') 225 RETURN 226 END 227C /* Deck rsplan */ 228 SUBROUTINE RSPLAN(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 229C 230C Lagrangian solution vector is returned in WRK(1) 231C 232#include "implicit.h" 233#include "dummy.h" 234#include "iratdef.h" 235#include "inftap.h" 236C 237 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 238 DIMENSION XINDX(*),WRK(*) 239C 240#include "priunit.h" 241#include "pgroup.h" 242#include "infpri.h" 243#include "infrsp.h" 244#include "wrkrsp.h" 245#include "rspprp.h" 246#include "inflr.h" 247#include "infdim.h" 248#include "inforb.h" 249#include "infesr.h" 250C 251C Local variables 252C 253 CHARACTER*8 LABEL, BLANK 254 PARAMETER ( D0 = 0.0D0 , DM8 = 1.0D-8, BLANK = ' ' ) 255C 256 CALL QENTER('RSPLAN') 257C 258C DETERMINE SECOND ORDER MOLECULAR PROPERTIES 259C 260 KREDE = 1 261 KREDS = KREDE + MAXRM*MAXRM 262 KIBTYP = KREDS + MAXRM*MAXRM 263 KEIVAL = KIBTYP + MAXRM 264 KRESID = KEIVAL + MAXRM 265 KEIVEC = KRESID + MAXRM 266 KREDGP = KEIVEC + MAXRM*MAXRM 267 KGP = KREDGP + MAXRM 268 KWRK1 = KGP + KZYVAR 269 LWRK1 = LWRK - KWRK1 270C 271 IF (LWRK1.LT.3*KZYVAR) THEN 272 WRITE (LUERR,9100) LWRK1,3*KZYVAR 273 CALL QTRACE(LUERR) 274 CALL QUIT('RSPLAN ERROR, INSUFFICIENT MEMORY') 275 ENDIF 276C 277C WORK SPACE FOR RSPEVE 278C 279 KWRKE = KWRK1 280 KBVECS = KWRKE + KZYVAR 281 9100 FORMAT(/' RSPLAN, work space too small for 3 (Z,Y)-vectors', 282 * /' had:',I10,', need more than:',I10) 283C 284 KZRED = 0 285 KZYRED = 0 286 THCRSP = THCESR 287 IPRRSP = IPRESR 288 MAXIT = MAXESR 289C 290C Call RSPCTL to solve linear set of response equations 291C 292 IF (TRPLET) THEN 293 LABEL = 'TRIPLAGR' 294 ELSE 295 LABEL = 'SINGLAGR' 296 END IF 297 WRITE (LUPRI,'(//A,I2,3A/2A)') 298 & ' RSPLAN -- linear response calculation for symmetry',KSYMOP, 299 & ' ( ',REP(KSYMOP-1),')', 300 & ' RSPLAN -- operator label : ',LABEL 301 CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM,WRK(KGP),LWRK1) 302 IF (IPRRSP.GT.120) THEN 303 WRITE(LUPRI,'(/A)') ' RSPLAN: LAGRANGIAN VECTOR ' 304 CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 305 END IF 306 XNOR = DDOT(KZYVAR,WRK(KGP),1,WRK(KGP),1) 307 IF ( XNOR.LT.DM8 ) THEN 308 WRITE(LUPRI,'(/A)') 309 * ' LAGRANGIAN VECTOR IS ZERO, SOLUTION IS SET TO ZERO' 310 CALL DZERO(WRK(1),KZYVAR) 311 GO TO 9999 312 END IF 313 WRK(KEIVAL) = D0 314 KEXSIM = 1 315 KEXCNV = 1 316 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 317 * .TRUE.,LABEL,BLANK,WRK(KGP),WRK(KREDGP), 318 * WRK(KREDE),WRK(KREDS), 319 * WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC), 320 * XINDX,WRK(KWRK1),LWRK1) 321C CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, 322C * LINEQ,GP,REDGP,REDE,REDS, 323C * IBTYP,EIVAL,EIVEC,XINDX,WRK,LWRK) 324C 325 NBX = 1 326 IBOFF = 0 327 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), 328 * WRK(KBVECS),WRK(KWRKE),NBX,IBOFF) 329C CALL RSPEVE(IBTYP,EIVAL,EIVEC,BVECS,WRK,NBX,IBOFF) 330 IF ( IPRRSP.GE.10 ) THEN 331 WRITE (LUPRI,'(/A)') ' SOLUTION VECTOR FOR LAGRANGIAN ' 332 CALL RSPPRC(WRK(KBVECS),KZCONF,KZVAR,LUPRI) 333 CALL RSPPRO(WRK(KBVECS+KZCONF),KZVAR,UDV,LUPRI) 334 END IF 335 CALL DCOPY(KZYVAR,WRK(KBVECS),1,WRK(1),1) 336 IF ( IPRRSP.GE.10 ) THEN 337 CALL HEADER(LABEL, LUPRI) 338 CALL OUTPUT(WRK, 1, KZVAR, 1, 2, KZVAR, 2, 1, LUPRI) 339 END IF 340 CALL WRTRSP( & 341 & LURSP, KZYVAR, WRK, LABEL, ' ', D0, D0, 1, 1, D0, .TRUE.& 342 &) 343C 344C *** end of RSPLAN -- 345C 346 9999 CALL QEXIT('RSPLAN') 347 RETURN 348 END 349C /* Deck esrinp */ 350 SUBROUTINE ESRINP(WORD) 351C 352C Module for reading "*ESR " input under **RESPONSE 353C 354#include "implicit.h" 355C 356#include "priunit.h" 357#include "gnrinf.h" 358#include "infrsp.h" 359#include "rspprp.h" 360#include "inflr.h" 361#include "infesr.h" 362#include "infpri.h" 363#include "zfs.h" 364#include "gtensor.h" 365#include "maxorb.h" 366#include "infinp.h" 367#include "parnmr.h" 368C 369 LOGICAL NEWDEF 370 PARAMETER ( NTABLE = 10 ) 371 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 372 CHARACTER*8 LABEL 373 CHARACTER LABDAT(9*ATMNUM)*8 374C 375 DATA TABLE /'.TRPPRP', '.SNGPRP', '.MAX IT', '.THCESR', 376 * '.PRINT ', '.G-TENS', '.ZFS ', '.FCCALC', 377 * '.ATOMS ','.SDCALC' / 378C 379C Initialize new HFC input 380 FCFLG = .FALSE. 381 SDFLG = .FALSE. 382 ESRNUC = 0 383 DO 50 I = 1, ATMNUM 384 DO 55 IH = 1, 7 385 ISODAT(IH,I) = 0 386 55 CONTINUE 387 50 CONTINUE 388C 389 NEWDEF = (WORD .EQ. '*ESR ') 390 ICHANG = 0 391 ESRCAL = .FALSE. 392 GCALC = .FALSE. 393 ZFSCAL = .FALSE. 394 IF (NEWDEF) THEN 395 WORD1 = WORD 396 100 CONTINUE 397 READ (LUCMD, '(A7)') WORD 398 CALL UPCASE(WORD) 399 PROMPT = WORD(1:1) 400 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') GO TO 100 401 IF (PROMPT .EQ. '.') THEN 402 ICHANG = ICHANG + 1 403 DO 200 I = 1, NTABLE 404 IF (TABLE(I) .EQ. WORD) THEN 405 GO TO (1,2,3,4,5,6,7,8,9,10), I 406 END IF 407 200 CONTINUE 408 IF (WORD .EQ. '.OPTION') THEN 409 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 410 GO TO 100 411 END IF 412 WRITE (LUPRI,'(/,3A,/)') ' KEYWORD "',WORD, 413 * '" NOT RECOGNIZED IN ESRINP.' 414 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 415 CALL QUIT(' ILLEGAL KEYWORD IN LRINP ') 416 1 CONTINUE ! .TRPPRP 417 READ(LUCMD,'( A )')LABEL 418 ESROPT( INDPRP(LABEL)) = .TRUE. 419 GO TO 100 420 2 CONTINUE ! .SNGPRP 421 READ(LUCMD,'( A )')LABEL 422 ESROPS( INDPRP(LABEL)) = .TRUE. 423 GO TO 100 424 3 CONTINUE ! .MAX IT 425 READ (LUCMD,*) MAXESR 426 GO TO 100 427 4 CONTINUE ! .THCESR 428 READ (LUCMD,*) THCESR 429 GO TO 100 430 5 CONTINUE ! .PRINT 431 READ (LUCMD,*) IPRESR 432 GO TO 100 433 6 CONTINUE ! .G-TENSOR 434 CALL GINP(WORD) 435 GO TO 100 436 7 CONTINUE ! .ZFS 437 ZFSCAL = .TRUE. 438 GO TO 100 439 8 CONTINUE ! .FCCALC 440 FCFLG = .TRUE. 441 GO TO 100 442 9 CONTINUE ! .ATOMS 443 READ (LUCMD,*) ESRNUC 444 IF (ESRNUC .GT. ATMNUM) THEN 445 WRITE(LUPRI,'(/A,I4,A,I4/A)') 446 & ' Sorry, but .ATOMS under *ESR is',ESRNUC, 447 & ' which is bigger than ATMNUM in parnmr.h',ATMNUM, 448 & ' Either reduce .ATOMS or increase ATMNUM and recompile.' 449 CALL QUIT('.ATOMS too big, see output') 450 END IF 451 READ (LUCMD,*) (NUCINF(IG), IG = 1, ESRNUC) 452 GO TO 100 453 10 CONTINUE ! .SDCALC 454 SDFLG = .TRUE. 455 GO TO 100 456 ELSE IF (PROMPT .EQ. '*') THEN 457 GO TO 300 458 ELSE 459 WRITE (LUPRI,'(/3A/)') ' PROMPT "',WORD, 460 * '" NOT RECOGNIZED IN ESRINP FOR *ESR.' 461 CALL QUIT(' ILLEGAL PROMPT under *ESR') 462 END IF 463 GO TO 100 464 END IF 465 300 CONTINUE 466 IF (THR_REDFAC .GT. 0.0D0) THEN 467 ICHANG = ICHANG + 1 468 WRITE (LUPRI,'(3A,1P,D10.2)') '@ INFO ',WORD1, 469 & ' thresholds multiplied with general factor',THR_REDFAC 470 THCESR = THCESR*THR_REDFAC 471 END IF 472C 473C FC labels generation 474C 475 IF (FCFLG) THEN 476 DO I = 1, ESRNUC 477 CALL FCOPER(NUCINF(I),LABEL) 478 ESROPT( INDPRP(LABEL)) = .TRUE. 479 END DO 480 IF (.NOT. SDFLG) CALL SETDEG 481 END IF 482C 483C SD labels generation 484C 485 IF (SDFLG) THEN 486 DO I = 1, ESRNUC 487 CALL SDOPER(NUCINF(I),LABDAT,INUM) 488 DO IG = 1, INUM 489 ESROPT( INDPRP(LABDAT(IG))) = .TRUE. 490 END DO 491 END DO 492 END IF 493C 494C Check degen array 495C 496 DO I = 1, ESRNUC 497 IF (DEGEN(I) .LT. 1.0) DEGEN(I)=1.0 498 END DO 499C 500C Setting isotopes data for HFC/pNMR printing 501C 502 CALL SETISO 503C 504 WRITE(LUPRI,'(/A)')' ********* ESRINP ********' 505C 506C 507 NTOESR = 0 508 DO I = 1,NPRLBL 509 IF (ESROPS(I)) NTOESR = NTOESR + 1 510 IF (ESROPT(I)) NTOESR = NTOESR + 1 511 END DO 512 ESRCAL = NTOESR .GT. 0 .OR. FCFLG .OR. SDFLG 513 IF (ESRCAL .OR. GCALC .OR. ZFSCAL) THEN 514 CALL HEADER('Changes of defaults under *ESR :',0) 515 IF ( ESRCAL ) THEN 516 WRITE (LUPRI,'(/A)') 517 * ' ESR - hyperfine calculation carried out' 518 IF (MAXESR.NE.60) WRITE(LUPRI,'(/A,I6)') 519 * ' MAXIMUM NUMBER OF ITERATIONS. MAXESR =',MAXESR 520 IF (THCESR.NE.1.0D-5) WRITE(LUPRI,'(/A,D13.6)') 521 * ' THRESHOLD FOR CONVERGENCE. THCESR =',THCESR 522 IF (IPRESR.NE.2) WRITE(LUPRI,'(/A,I5)') 523 * ' PRINT LEVEL IN ESR ROUTINES. IPRESR =',IPRESR 524 END IF 525 IF (ZFSCAL) THEN 526 WRITE (LUPRI,'(/A)') 527 * ' ESR - zero-field splitting calculation carried out' 528 END IF 529 IF (GCALC) THEN 530 WRITE (LUPRI,'(/A)') 531 * ' ESR - g-tensor calculation carried out' 532 IF (ECC) WRITE(LUPRI,'(A)') 533 * ' Electronic charge centroid used as gauge origin' 534 IF (MNFSO) WRITE(LUPRI,'(A)') 535 * ' Mean-field spin-orbit approximation is used' 536 IF (SCALED_CHARGES) WRITE(LUPRI,'(A)') 537 * ' Two-electron spin-orbit by scaled nuclear charges' 538 IF (.NOT.DOALL) THEN 539 WRITE(LUPRI,'(/A)') 540 * ' Selected contributions to electronic g_tensor:' 541 IF (DORMC) WRITE(LUPRI,'(/A)') 542 * ' Relativistic mass correction' 543 IF (DOGC1) WRITE(LUPRI,'(/A)') 544 * ' One-electron gauge correction' 545 IF (DOGC2) WRITE(LUPRI,'(/A)') 546 * ' Two-electron gauge correction' 547 IF (DOOZSO1) WRITE(LUPRI,'(/A)') 548 * ' One-electron spin-orbit + orbital Zeeman correction' 549 IF (DOOZSO2) WRITE(LUPRI,'(/A)') 550 * ' Two-electron spin-orbit + orbital Zeeman correction' 551 END IF 552 END IF 553 ELSE 554 WRITE (LUPRI,'(//A/)') 555 * 'INFO: "*ESR " input ignored because no operators requested.' 556 END IF 557 IF (ISPIN .EQ. 1) THEN 558 NWARN = NWARN + 1 559 WRITE (LUPRI,'(//A/)') 560 * ' WARNING: "*ESR " input ignored because'// 561 & ' singlet reference state has no ESR properties !!!' 562 ESRCAL = .FALSE. 563 FCFLG = .FALSE. 564 SDFLG = .FALSE. 565 GCALC = .FALSE. 566 ZFSCAL = .FALSE. 567 END IF 568C 569C *** END OF ESRINP 570C 571 RETURN 572 END 573C /* Deck rspesr */ 574 SUBROUTINE RSPESR(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 575C 576C Calculate triplet expectation values with Lagrangian correction 577C and singlet expectation values for open shell systems with MCSCF or CI 578C (non-zero spin densities). 579C 580C Revised March 2003 hjaaj 581C 582#include "implicit.h" 583#include "iratdef.h" 584C 585 DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) 586 DIMENSION XINDX(*),WRK(*) 587C 588 PARAMETER ( D0 = 0.0D0, DUMMY = 1.0D+20 ) 589 LOGICAL TRPSAVE 590C 591#include "maxorb.h" 592#include "infinp.h" 593C 594#include "priunit.h" 595#include "infpri.h" 596#include "infrsp.h" 597#include "wrkrsp.h" 598#include "rspprp.h" 599#include "inflr.h" 600#include "infdim.h" 601#include "inforb.h" 602#include "infesr.h" 603C 604C 605C Common blocks for HFC/pNMR printing 606C 607#include "codata.h" 608#include "gfac.h" 609#include "mxcent.h" 610#include "nuclei.h" 611#include "chrxyz.h" 612#include "parnmr.h" 613C 614C Local variables for HFC values handling 615C 616 CHARACTER*8 LABEL 617 CHARACTER*2 CTMP 618C 619 CALL QENTER('RSPESR') 620 CALL HEADER('Output from RSPESR module',0) 621C HFC printing counters 622 IFCIND = 0 623 REFSPIN=DBLE(ISPIN-1) 624C 625C Zero SDVAL 626C 627 DO 50 I = 1, ATMNUM 628 DO 55 IH = 1, 3 629 DO 60 IG = 1, 3 630 SDVAL(IG,IH,I) = 0.0D0 631 60 CONTINUE 632 55 CONTINUE 633 50 CONTINUE 634C 635C DETERMINE SECOND ORDER MOLECULAR PROPERTIES 636C 637 IF ( NESRT(KSYMOP).GT. 0 ) THEN 638 NDREF = KZCONF 639C ... we use determinants when triplet, thus not NCREF 640 NDREF = MAX(KZCONF,NCREF) 641Chj ... we use determinants triplet, thus KZCONF, 642Chj ... However, for ROHF NCREF = 1, KZCONF = 0 and 643Chj we need CREF(1) for RSPDM. Thus MAX(KZCONF,NCREF) 644C 645C Note that GETGPV uses WRK(KGP) as scratch, thus KGP last allocation ! 646C 647 IF (RSPCI) THEN 648 KTUDV = 1 649 KCREF = KTUDV + N2ASHX 650 KWRK1 = KCREF + NDREF 651 KLAGR = KCREF 652 KGP = KCREF 653 ELSE 654 KLAGR = 1 655 CALL RSPLAN(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) 656 IF (IPRRSP.GT.120) THEN 657 WRITE(LUPRI,'(/A)') 658 & 'RSPESR: solution vector for lagrangian :' 659 CALL OUTPUT(WRK,1,KZVAR,1,2,KZVAR,2,1,LUPRI) 660 END IF 661 KTUDV = KLAGR + KZYVAR 662 KCREF = KTUDV + N2ASHX 663 KGP = KCREF 664 KWRK1 = KGP + MAX(KZYVAR,NDREF) 665 END IF 666 LWRK1 = LWRK - KWRK1 667 IF (LWRK1.LT.0) CALL ERRWRK('RSPESR',KWRK1-1,LWRK) 668C 669C Get one electron spin density (triplet MS=0 density) 670C 671 CALL GETREF(WRK(KCREF),NDREF) 672 KFREE = 1 673 LFREE = LWRK1 674 ISPIN1 = 1 675 ISPIN2 = 0 676 CALL RSPDM(IREFSY,IREFSY,NDREF,NDREF,WRK(KCREF),WRK(KCREF), 677 * WRK(KTUDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 678 * XINDX,WRK(KWRK1),KFREE,LFREE) 679C CALL RSPDM(ISYM,ISYM,NCDIM.NCDIM,REF,REF, 680C TUDV,DUMMY, ISPIN1,ISPIN2,TDM,NORHO2, 681C XINDX,WRK,KFREE,LFREE) 682 IF (IPRRSP.GE.5) THEN 683 WRITE(LUPRI,'(/A)')'RSPESR: one electron spin density ' 684 CALL OUTPUT(WRK(KTUDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 685 END IF 686 TRPSAVE= TRPLET 687 TRPLET = .TRUE. 688 DO 100 IOP = 1,NESRT(KSYMOP) 689 IF (.NOT. RSPCI ) THEN 690 IF (IPRRSP.GT.150) THEN 691 WRITE(LUPRI,'(/A)') 692 * 'RSPESR: Lagrangian vector before product' 693 CALL OUTPUT(WRK(KLAGR),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 694 END IF 695 CALL GETGPV(LBESRT(KSYMOP,IOP),FC,FV,CMO,UDV,PV,XINDX, 696 * ANTSYM,WRK(KGP),LWRK1 ) 697C CALL GETGPV(LABELOP,FC,CMO,UDV,PV,XINDX,ANTSYM,WRK,LWRK) 698 IF (IPRRSP.GT.120) THEN 699 WRITE(LUPRI,'(/2A)') 700 * 'RSPESR: GP vector with label: ',LBESRT(KSYMOP,IOP) 701 CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) 702 END IF 703 FOPLG = -DDOT(KZYVAR,WRK(1),1,WRK(KGP),1) 704 IF (IPRRSP.GT.5) WRITE(LUPRI,'(/A,1P,D18.10)') 705 & ' GP * LAGRANG SOLUTION:',FOPLG 706C 707 ELSE 708 FOPLG = D0 709 END IF 710 CALL PRP1AVE(LBESRT(KSYMOP,IOP),AVEVAL,CMO, 711 * WRK(KTUDV),WRK(KWRK1),LWRK1,IPRRSP) 712 HYPFIN = AVEVAL + FOPLG 713 WRITE(LUPRI,'(/2A,3(A,F10.6))') 714 * 'TRIPLET OPERATOR: "',LBESRT(KSYMOP,IOP), 715 * '" LAGRANGIAN:',FOPLG,' AVERAGE:',AVEVAL,' TOTAL:',HYPFIN 716C 717C FC and SD values arrays formation 718C 719 IF (LBESRT(KSYMOP,IOP)(1:3) .EQ. 'FC ') THEN 720 IFCIND=IFCIND+1 721 HYPVAL(IFCIND,1)=AVEVAL 722 HYPVAL(IFCIND,2)=FOPLG 723 END IF 724 IF (LBESRT(KSYMOP,IOP)(1:3) .EQ. 'SD ') THEN 725 LABEL=LBESRT(KSYMOP,IOP) 726 READ (LABEL,'(3X,I3)') INDSD1 727 CTMP=LABEL(7:8) 728 indsd2 = -9898 729 IF (CTMP .EQ. ' x') INDSD2=1 730 IF (CTMP .EQ. ' y') INDSD2=2 731 IF (CTMP .EQ. ' z') INDSD2=3 732 if (indsd2 .eq. -9898) then 733 write(lupri,*) 'hjaaj label=',label, indsd1 734 write(lupri,*) 'hjaaj no x or y or z in label(7:8) !' 735 call quit('hjaaj: problem in lagrang.F') 736 end if 737 SDVAL(SDIND(2,INDSD1),INDSD2,SDIND(1,INDSD1))=HYPFIN 738 END IF 739C 740100 CONTINUE 741 TRPLET = TRPSAVE 742 END IF 743 DO 300 ISYM = 2,NSYM 744 DO 350 IOP = 1,NESRT(ISYM) 745 WRITE(LUPRI,'(/3A/A,I3)') 'TRIPLET OPERATOR: "', 746 & LBESRT(ISYM,IOP),'" contribution = 0.0 by symmetry.', 747 & '- Operator is of symmetry no.',ISYM 748350 CONTINUE 749300 CONTINUE 750C 751C Singlet operators: 752C 753 TRPSAVE= TRPLET 754 TRPLET = .FALSE. 755C ... so PRP1AVE calculates singlet expectation values /hjaaj march 2003 ... 756C (if TRPLET true, then inactive density matrix is omitted!) 757 DO 400 IOP = 1,NESRS(KSYMOP) 758 CALL PRP1AVE(LBESRS(KSYMOP,IOP),AVEVAL,CMO, 759 * UDV,WRK(KWRK1),LWRK1,IPRRSP) 760C CALL PRP1AVE(LABELOP,AVEVAL,CMO, UDV,WRK,LWRK,IPRINT) 761 WRITE(LUPRI,'(/3A,F10.6)') 762 * 'SINGLET OPERATOR: "',LBESRS(KSYMOP,IOP),'" AVERAGE:',AVEVAL 763400 CONTINUE 764C 765 DO 500 ISYM = 2,NSYM 766 DO 550 IOP = 1,NESRS(ISYM) 767 WRITE(LUPRI,'(/3A/A,I3)') 'SINGLET OPERATOR: "', 768 & LBESRS(ISYM,IOP),'" contribution = 0.0 by symmetry.', 769 & '- Operator is of symmetry no.',ISYM 770550 CONTINUE 771500 CONTINUE 772 TRPLET = TRPSAVE 773C 774C Isotropic hyperfine coupling printing 775C 776 IF (FCFLG) THEN 777 CALL AROUND('Isotropic Hyperfine Coupling') 778 WRITE(LUPRI,'(A)') ' Atom Mass G-val ' 779 * //' A, Mhz A, G ' 780 WRITE(LUPRI,'(A)') ' -----------------------------' 781 * //'--------------------------' 782 DO 600 ITMP = 1, ESRNUC 783 DO 610 ISONM = 1, ISODAT(1,NUCINF(ITMP)) 784 IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) 785 XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'GVAL') 786 XISMAS=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'MASS') 787 HYPFIN=((HYPVAL(ITMP,1)+HYPVAL(ITMP,2)) 788 * *XATGV*XTHZ*1.0D-6*ALPHA2) 789 * /(2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) 790 WRITE(LUPRI,'(A,F6.2,A,F8.5,A,F9.4,A,F9.4,A)') 791 * ' * '//NAMN(NUCINF(ITMP))(1:4)//' * ', XISMAS, 792 * ' * ', XATGV,' * ',HYPFIN,' * ',HYPFIN/2.8025D0, 793 * ' *' 794 610 CONTINUE 795 600 CONTINUE 796 WRITE(LUPRI,'(A)') ' -----------------------------' 797 * //'--------------------------' 798 WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' 799 * //' symmetry unique atoms!' 800 CALL AROUND('R-U contributions to Isotropic hyperfine' 801 * //' coupling') 802 WRITE(LUPRI,'(A)') ' Atom G-val Average, MHz ' 803 * //' Response, MHz Total, MHz ' 804 WRITE(LUPRI,'(A)') ' ------------------------------' 805 * //'---------------------------------' 806 DO 601 ITMP = 1, ESRNUC 807 DO 611 ISONM = 1, ISODAT(1,NUCINF(ITMP)) 808 IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) 809 XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'GVAL') 810 HYPFIN=(XATGV*XTHZ*1.0D-6*ALPHA2) 811 * /(2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) 812 WRITE(LUPRI,'(A,F8.5,A,F11.4,A,F11.4,A,F11.4,A)') 813 * ' * '//NAMN(NUCINF(ITMP))(1:4)//' * ', XATGV, 814 * ' * ',HYPFIN*HYPVAL(ITMP,1),' * ', 815 * HYPFIN*HYPVAL(ITMP,2),' * ', 816 * HYPFIN*(HYPVAL(ITMP,1)+HYPVAL(ITMP,2)),' * ' 817 611 CONTINUE 818 601 CONTINUE 819 WRITE(LUPRI,'(A)') ' ------------------------------' 820 * //'---------------------------------' 821 WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' 822 * //' symmetry unique atoms!' 823 ENDIF 824C 825C Anisotropic hyperfine coupling printing 826C 827 IF (SDFLG) THEN 828 CALL AROUND('Anisotropic Hyperfine Coupling') 829 DO 700 ITMP=1,ESRNUC 830 DO 710 ISONM=1,ISODAT(1,NUCINF(ITMP)) 831 IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) 832 XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS, 833 * 'GVAL') 834 XISMAS=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS, 835 * 'MASS') 836 HYPFIN=(XATGV*XTHZ*1.0D-6*ALPHA2)/ 837 * (2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) 838 WRITE(LUPRI,'(/A,F6.2,A,F8.5,A)') 839 * ' A Tensor components for '// 840 * NAMN(NUCINF(ITMP))(1:4)//' ( Mass =', 841 * XISMAS, ' G-val =', XATGV,') in MHz: ' 842 WRITE(LUPRI,'(A)') ' =================================' 843 * //'==================================' 844 WRITE(LUPRI,'(/A)') ' ' 845 WRITE(LUPRI,'(A)') ' SX SY SZ ' 846 * //' ' 847 DO 720 ICRD=1,3 848 WRITE(LUPRI,'(A,F11.4,A,F11.4,A,F11.4,A)') 849 * ' I'//CHRXYZ(ICRD)//' * ', 850 * HYPFIN*SDVAL(ICRD,1,NUCINF(ITMP)), ' * ', 851 * HYPFIN*SDVAL(ICRD,2,NUCINF(ITMP)), ' * ', 852 * HYPFIN*SDVAL(ICRD,3,NUCINF(ITMP)), ' * ' 853 720 CONTINUE 854 710 CONTINUE 855 700 CONTINUE 856 WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' 857 * //' symmetry unique atoms!' 858 ENDIF 859C 860C *** end of RSPESR -- 861 CALL QEXIT('RSPESR') 862 RETURN 863 END 864C /* Deck solgdt */ 865 SUBROUTINE SOLGDT(CREF,CMO,INDXCI,DV,G,ESOLT,ERLM,WRK, 866 & LFREE,NHCREF,KREFSY) 867C 868C Based on SOLGRD: 869C Copyright 29-Nov-1986 Hans Joergen Aa. Jensen 870C 871C Purpose: calculate MCSCF energy and gradient contribution 872C from a surrounding medium, cavity radius = Rsol 873C and dielectric constant = EPsol. 874C 875C Output: 876C G MCSCF gradient with solvation contribution added 877C ESOLT total solvation energy 878C ERLM(lm,1) contains Esol(l,m) contribution to ESOLT 879C ERLM(lm,2) contains Tsol(l,m) 880C 881#include "implicit.h" 882#include "dummy.h" 883C 884 DIMENSION CREF(*), CMO(*), INDXCI(*) 885 DIMENSION DV(*), G(*), ERLM(NLMSOL,2), WRK(*) 886 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0 ) 887#include "thrzer.h" 888#include "iratdef.h" 889#include "priunit.h" 890#include "infrsp.h" 891C 892C Used from common blocks: 893C INFINP: NLMSOL, LSOLMX, EPSOL, RSOL(3) 894C INFVAR: NCONF, NWOPT, NVAR, NVARH 895C INFORB: NNASHX, NNBASX, NNORBX, etc. 896C INFIND: IROW(*) 897C INFTAP: LUSOL, LUIT2 898C INFPRI: IPRSOL 899C 900#include "maxash.h" 901#include "maxorb.h" 902#include "infinp.h" 903#include "infvar.h" 904#include "inforb.h" 905#include "infind.h" 906#include "inftap.h" 907#include "infpri.h" 908C 909 LOGICAL FIRST 910 PARAMETER (MXLMAX = 50) 911 DIMENSION ISYRLM(2*MXLMAX+1) 912 CHARACTER*8 STAR8, SOLVDI, EODATA 913 SAVE FIRST 914 DATA FIRST/.TRUE./, STAR8/'********'/ 915 DATA SOLVDI/'SOLVDIAG'/, EODATA/'EODATA '/ 916C 917C Statement functions; 918C define automatic arrays (dynamic core allocation) 919C 920 FLVEC(LM) = WRK(LM) 921 FLINR(LM) = WRK(KFLINR-1+LM) 922 TLMSI(LM) = WRK(KTLMSI-1+LM) 923C 924 CALL QENTER('SOLGDT') 925C 926 IF (LSOLMX .GT. MXLMAX) THEN 927 WRITE (LUERR,*) 'ERROR SOLGDT, increase MXLMAX parameter' 928 WRITE (LUERR,*) ' LSOLMX =',LSOLMX 929 WRITE (LUERR,*) ' MXLMAX =',MXLMAX 930 CALL QUIT('ERROR SOLGDT, increase MXLMAX parameter') 931 END IF 932C 933C Core allocation 934C FLVEC f(l) factors in solvent energy expression 935C DIASH diagonal of solvent contribution to Hessian 936C GRDLM TELM gradient for current l,m value in the l,m loop 937C UCMO CMO unpacked (i.e. no symmetry blocking) 938C RLMAC active-active subblock of RLM 939C RLM R(l,m) integrals for current l,m value in l,m loop 940C 941C If (INERSF) 942C (i.e. If (inertial polarization contribution to final state)) 943C FLINR f(l) factors for inertial pol. contrib. 944C TLMSI T(lm) values for initial state 945C end if 946C 947 KFLVEC = 1 948C ... NOTE: KFLVEC = 1 assumed in FLVEC(LM) definition above. 949 IF (INERSF) THEN 950 KFLINR = KFLVEC + NLMSOL 951 KTLMSI = KFLINR + NLMSOL 952 KDIASH = KTLMSI + NLMSOL 953 ELSE 954 KFLINR = 1 955 KTLMSI = 1 956 KDIASH = KFLVEC + NLMSOL 957 END IF 958 KGRDLM = KDIASH + NVAR 959 KUCMO = KGRDLM + NVARH 960 KRLMAC = KUCMO + NORBT*NBAST 961 KRLM = KRLMAC + NNASHX 962 KW10 = KRLM + NNORBX 963C 1.1 read rlmao in ao basis and transform to rlm in mo basis 964 KRLMAO = KW10 965 KW20 = KRLMAO + NNBASX 966C 1.2 diagonal contribution for current l,m value in the l,m loop 967 KDIALM = KW10 968 KW21 = KDIALM + NVAR 969 LW21 = LFREE - KW21 970C 1.3 rest of CSF contribution 971 KW22 = KW10 972C 973 KTDV = MAX(KW20,KW21,KW22) 974 KWRK1 = KTDV + NASHT * NASHT 975 LWRK1 = LFREE - KWRK1 976 IF (LWRK1 .LT. 0) CALL ERRWRK('SOLGDT',-KWRK1,LFREE) 977C 978 IF (IPRSOL .GE. 130) THEN 979 WRITE (LUPRI,'(/A/A,2I10)') 980 * ' SOLGDT - gtot (input) - non-zero elements', 981 * ' NCONF, NWOPT =',NCONF,NWOPT 982 DO 40 I = 1,NCONF 983 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 984 * ' conf #',I,G(I) 985 40 CONTINUE 986 DO 50 I = NCONF+1,NVAR 987 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 988 * ' orb #',I,G(I) 989 50 CONTINUE 990 END IF 991C 992C Calculate f(l) factors 993C If (INERSF) FLVEC factors describe the optical polarization 994C and FLINR factors describe the inertial polarization 995C else FLVEC may describe optical or static polarization 996C 997 CALL SOLFL(WRK(KFLVEC),EPSOL,RSOL,LSOLMX) 998 IF (INERSF) THEN 999 CALL SOLINR(WRK(KFLVEC),WRK(KFLINR),WRK(KTLMSI)) 1000 END IF 1001 IF ((IPRSOL .GE. 5 .AND. FIRST) .OR. IPRSOL .GE. 15) THEN 1002 IF (.NOT.INERSF) THEN 1003 WRITE (LUPRI,'(//A/A)') 1004 * ' SOLGDT: l f(l) factor', 1005 * ' === =================' 1006 ELSE 1007 WRITE (LUPRI,'(//A/A)') 1008 * ' SOLGDT: l optical f(l) factor inertial f(l) factor', 1009 * ' === =================== ====================' 1010 END IF 1011 DO 140 L = 0,LSOLMX 1012 LL = (L+1)*(L+1) 1013 FL = FLVEC(LL) 1014 IF (INERSF) THEN 1015 FLI = FLINR(LL) 1016 WRITE (LUPRI,'(I15,F17.10,F21.10)') L, FL, FLI 1017 ELSE 1018 WRITE (LUPRI,'(I15,F16.10)') L, FL 1019 END IF 1020 140 CONTINUE 1021 END IF 1022C 1023C Read and check dimension information (if first read) and 1024C nuclear contributions to ERLM (always). 1025C 1026 CALL GPOPEN(LUSOL,FNSOL,'OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) 1027 REWIND LUSOL 1028 CALL MOLLAB('SOLVRLM ',LUSOL,LUERR) 1029 IF (FIRST) THEN 1030 READ (LUSOL) LMAXSS, LMTOT, NNNBAS 1031 NERR = 0 1032 IF (LMAXSS .LT. LSOLMX) THEN 1033 NERR = NERR + 1 1034 WRITE (LUPRI,'(//2A,2(/A,I5))') ' SOLGDT ERROR,', 1035 * ' insufficient number of integrals on LUSOL', 1036 * ' l max from SIRIUS input :',LSOLMX, 1037 * ' l max from LUSOL file :',LMAXSS 1038 END IF 1039 IF ((LMAXSS+1)**2 .NE. LMTOT) THEN 1040 NERR = NERR + 1 1041 WRITE (LUPRI,'(//2A,3(/A,I5))') ' SOLGDT ERROR,', 1042 * ' LUSOL file info inconsistent', 1043 * ' l_max :',LMAXSS, 1044 * ' (l_max + 1) ** 2 :',(LMAXSS+1)**2, 1045 * ' LMTOT :',LMTOT 1046 END IF 1047 IF (NNNBAS .NE. NBAST) THEN 1048 NERR = NERR + 1 1049 WRITE (LUPRI,'(//2A,3(/A,I5))') ' SOLGDT ERROR,', 1050 * ' LUSOL file info inconsistent with SIRIUS input', 1051 * ' NBAST - LUSOL :',NNNBAS, 1052 * ' NBAST - SIRIUS :',NBAST 1053 END IF 1054 IF (NERR .GT. 0) THEN 1055 CALL QUIT('SOLGDT ERROR: LUSOL file not OK for this calc.') 1056 END IF 1057 ELSE 1058 READ (LUSOL) 1059 END IF 1060 CALL READT(LUSOL,NLMSOL,ERLM(1,2)) 1061C 1062 IF (IPRSOL .GE. 20 .AND. NASHT .GT. 0) THEN 1063 WRITE (LUPRI,'(/A)') ' SOLGDT - DV matrix :' 1064 CALL OUTPAK(DV,NASHT,1,LUPRI) 1065 END IF 1066 IF (IPRSOL .GE. 7) THEN 1067 WRITE (LUPRI,'(/A/)') 1068 * ' l, m, Tn(lm) - the nuclear contributions :' 1069 LM = 0 1070 DO 220 L = 0,LSOLMX 1071 DO 210 M = -L,L 1072 LM = LM + 1 1073 WRITE (LUPRI,'(2I5,F15.10)') L,M,ERLM(LM,2) 1074 210 CONTINUE 1075 WRITE (LUPRI,'()') 1076 220 CONTINUE 1077 END IF 1078C 1079C Unpack symmetry blocked CMO 1080C Loop over l,m expansion 1081C 1082 CALL UPKCMO(CMO,WRK(KUCMO)) 1083 IF (IPRSOL .GE. 6) 1084 * WRITE (LUPRI, '(//A/)') ' SOLGDT: START LOOP OVER LM' 1085 CALL DZERO(WRK(KDIASH),NVAR) 1086 LM = 0 1087 DO 520 L = 0,LSOLMX 1088 READ (LUSOL) L1,(ISYRLM(M),M=1,2*L+1) 1089 IF (L1 .NE. L) THEN 1090 WRITE (LUERR,*) 'ERROR SOLGDT: L from LUSOL not as expected' 1091 WRITE (LUERR,*) 'L from 520 loop:',L 1092 WRITE (LUERR,*) 'L from LUSOL :',L1 1093 CALL QUIT('ERROR SOLGDT: L from LUSOL not as expected') 1094 END IF 1095 DO 500 M = -L,L 1096 LM = LM + 1 1097 IF (IPRSOL .GE. 15) THEN 1098 WRITE (LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, 1099 * ' ====================' 1100 WRITE (LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYRLM(L+M+1)) 1101 END IF 1102 IF (ISYRLM(L+M+1) .NE. 1) THEN 1103 IF (ABS(ERLM(LM,2)) .GT. 1000.D0*THRZER) THEN 1104 WRITE (LUPRI,*) 'ERROR SOLGDT for l,m',L,M 1105 WRITE (LUPRI,*) 'Symmetry :',ISYRLM(L+M+1) 1106 WRITE (LUPRI,*) 'Tn(l,m) .ne. 0, but =',ERLM(LM,2) 1107 CALL QUIT('ERROR SOLGDT: Tn(l,m) not 0 as expected') 1108 END IF 1109 ERLM(LM,2) = D0 1110C ... to fix round-off errors in Tn(l,m) calculation 1111 IF (ISYRLM(L+M+1) .GT. 1) READ (LUSOL) 1112 GO TO 500 1113 END IF 1114C 1115C Read R(l,m) in ao basis and transform to mo basis. 1116C Extract active-active block in RLMAC(1) = WRK(KRLMAC). 1117C 1118 CALL READT(LUSOL,NNBASX,WRK(KRLMAO)) 1119 CALL UTHU(WRK(KRLMAO),WRK(KRLM),WRK(KUCMO),WRK(KWRK1), 1120 & NBAST,NORBT) 1121 IF (NASHT .GT. 0) THEN 1122 CALL GETAC2(WRK(KRLM),WRK(KRLMAC)) 1123 END IF 1124 IF (IPRSOL .GE. 15) THEN 1125 WRITE (LUPRI,'(/A)') ' Rlm_ao matrix:' 1126 CALL OUTPAK(WRK(KRLMAO),NBAST,1,LUPRI) 1127 WRITE (LUPRI,'(/A)') ' Rlm_mo matrix:' 1128 CALL OUTPAK(WRK(KRLM), NORBT,1,LUPRI) 1129 IF (NASHT .GT. 0) THEN 1130 WRITE (LUPRI,'(/A)') ' Rlm_ac matrix:' 1131 CALL OUTPAK(WRK(KRLMAC),NASHT,1,LUPRI) 1132 END IF 1133 END IF 1134C 1135C Add electronic contribution TE(l,m) to T(l,m) 1136C 1137 KFREE=1 1138 ISPIN1=0 1139 ISPIN2=0 1140 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, 1141 * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 1142 * INDXCI,WRK(KWRK1),KFREE,LWRK1) 1143C 1144C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX 1145C 1146 CALL DSITSP(NASHT,WRK(KTDV),DV) 1147C 1148 TELM = SOLELM(DV,WRK(KRLMAC),WRK(KRLM),TELMAC) 1149C 1150C construct again triplet density matrix 1151 KFREE =1 1152 ISPIN1=1 1153 ISPIN2=0 1154 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, 1155 * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 1156 * INDXCI,WRK(KWRK1),KFREE,LWRK1) 1157C 1158C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX 1159C 1160 CALL DSITSP(NASHT,WRK(KTDV),DV) 1161C 1162 IF (IPRSOL .GE. 6) THEN 1163 WRITE (LUPRI,'(A,2I5,/A,3F17.8)') 1164 * ' --- l, m :',L,M, 1165 * ' Te(lm), Tn(lm), T(lm) :', 1166 * TELM,ERLM(LM,2),ERLM(LM,2)-TELM 1167 IF (IPRSOL .GE. 10) WRITE (LUPRI,'(A,F17.8)') 1168 * ' --- active part of Te(lm) :',TELMAC 1169 IF (INERSF) WRITE (LUPRI,'(A,F17.8)') 1170 * ' --- inertial T(lm) value :',TLMSI(LM) 1171 END IF 1172 ERLM(LM,2) = ERLM(LM,2) - TELM 1173 IF (ABS(ERLM(LM,2)) .LE. THRZER) THEN 1174 ERLM(LM,2) = D0 1175 GO TO 500 1176 END IF 1177C ... test introduced 880109 hjaaj 1178C (the only possible problem is the DO 420 loop, 1179C but I think (w.o. having checked) that this 1180C contribution to the Hessian diagonal also will be 1181C zero if ERLM(LM,2) zero). 1182C 1183C Calculate orbital TE(l,m) gradient contribution 1184C and part of csf contribution. 1185C 1186 CALL DZERO(WRK(KGRDLM),NVARH) 1187 IF (NCONF .GT. 1) THEN 1188 CALL SOLGC(CREF,WRK(KRLMAC),TELMAC,WRK(KGRDLM),INDXCI, 1189 & WRK(KWRK1),LWRK1) 1190C CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) 1191 END IF 1192 IF (NWOPT .GT. 0) THEN 1193 CALL SOLGO(D0,DV,WRK(KRLM),WRK(KGRDLM+NCONF)) 1194 END IF 1195C 1196C 1197C Obtain DIALM = diagonal TE(l,m) Hessian 1198C = 2 ( <i|R(l,m)|i> - TE(l,m) ) 1199C Add the DIALM contribution and the GRDLM contribution 1200C to solvent Hessian diagonal. 1201C 1202Clf the diagonal hessian is not needed for triplet gradients 1203Clf CALL SOLDIA(TELMAC,WRK(KRLMAC),INDXCI, 1204Clf * WRK(KRLM),DV,WRK(KDIALM),WRK(KW21),LW21) 1205C CALL SOLDIA(TELM,RLMAC,INDXCI,RLM,DV,DIAG,WRK,LFREE) 1206C 1207 FAC1 = - D2 * FLVEC(LM) * ERLM(LM,2) 1208 IF (INERSF) THEN 1209 FAC1 = FAC1 - FLINR(LM) * D2 * TLMSI(LM) 1210 END IF 1211 FAC2 = D2 * FLVEC(LM) 1212 DO 420 I = 0,(NVAR-1) 1213 WRK(KDIASH+I) = WRK(KDIASH+I) 1214 * + FAC1 * WRK(KDIALM+I) 1215 * + FAC2 * WRK(KGRDLM+I) * WRK(KGRDLM+I) 1216 420 CONTINUE 1217C 1218C test orthogonality 1219C 1220 IF (IPRSOL .GE. 120) THEN 1221 WRITE (LUPRI,'(/A)')' SOLGDT - grdlm, dialm, diash, cref' 1222 DO 430 I = 1,NCONF 1223 WRITE (LUPRI,'(A,I10,4F12.6)') ' conf #',I, 1224 * WRK(KGRDLM-1+I),WRK(KDIALM-1+I),WRK(KDIASH-1+I),CREF(I) 1225 430 CONTINUE 1226 END IF 1227 TEST = DDOT(NCONF,CREF,1,WRK(KGRDLM),1) 1228 IF (ABS(TEST) .GT. 1.D-8) THEN 1229 NWARN = NWARN + 1 1230 WRITE (LUPRI,'(/A,I5,/A,1P,D12.4)') 1231 * ' SOLGDT WARNING, for LM =',LM, 1232 * ' <CREF | GRADlm > =',TEST 1233 END IF 1234C 1235C Add TE(l,m) gradient contribution to MCSCF gradient 1236C g = g - 2 f(l) * T(l,m) * (dTE(l,m)/d(lambda)) 1237C 1238 FAC = - D2 * FLVEC(LM) * ERLM(LM,2) 1239 IF (INERSF) THEN 1240 FAC = FAC - FLINR(LM) * D2 * TLMSI(LM) 1241 END IF 1242 CALL DAXPY(NVARH,FAC,WRK(KGRDLM),1,G,1) 1243 IF (IPRSOL .GE. 140) THEN 1244 WRITE (LUPRI,'(/A/A,2I10)') 1245 * ' SOLGDT - grdlm, gtot (accum) - non-zero grdlm', 1246 * ' NCONF, NWOPT =',NCONF,NWOPT 1247 DO 440 I = 1,NCONF 1248 IF (WRK(KGRDLM-1+I) .NE. D0) 1249 * WRITE (LUPRI,'(A,I10,3F15.10)') 1250 * ' conf #',I,FAC*WRK(KGRDLM-1+I),G(I) 1251 440 CONTINUE 1252 DO 450 I = NCONF+1,NVAR 1253 IF (WRK(KGRDLM-1+I) .NE. D0) 1254 * WRITE (LUPRI,'(A,I10,3F15.10)') 1255 * ' orb #',I,FAC*WRK(KGRDLM-1+I),G(I) 1256 450 CONTINUE 1257 END IF 1258C 1259 500 CONTINUE 1260 520 CONTINUE 1261C 1262 CALL GPCLOSE(LUSOL,'KEEP') 1263C 1264C 500 is end of (l,m) loop. 1265C 1266C 1267 IF (IPRSOL .GE. 130) THEN 1268 WRITE (LUPRI,'(/A/A,2I10)') 1269 * ' SOLGDT - gtot (output) - non-zero elements', 1270 * ' NCONF, NWOPT =',NCONF,NWOPT 1271 DO 840 I = 1,NCONF 1272 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1273 * ' conf #',I,G(I) 1274 840 CONTINUE 1275 DO 850 I = NCONF+1,NVAR 1276 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1277 * ' orb #',I,G(I) 1278 850 CONTINUE 1279 END IF 1280C 1281C 1282C Calculate ER(l,m) energy contributions and add them up 1283C 1284 ESOLT = D0 1285 DO 900 LM = 1,NLMSOL 1286 ERLM(LM,1) = FLVEC(LM) * ERLM(LM,2) * ERLM(LM,2) 1287 IF (INERSF) THEN 1288 ERLM(LM,1) = ERLM(LM,1) 1289 * + FLINR(LM) * ERLM(LM,2) * D2 * TLMSI(LM) 1290 * - FLINR(LM) * TLMSI(LM) * TLMSI(LM) 1291 END IF 1292 ESOLT = ESOLT + ERLM(LM,1) 1293 900 CONTINUE 1294C 1295 FIRST = .FALSE. 1296 CALL QEXIT('SOLGDT') 1297 RETURN 1298C end of solgdt. 1299 END 1300C /* Deck FCOPER */ 1301 SUBROUTINE FCOPER(ATMIND,LABINT) 1302#include "implicit.h" 1303#include "priunit.h" 1304#include "mxcent.h" 1305 CHARACTER*8 LABINT 1306 INTEGER ATMIND 1307#include "nuclei.h" 1308#include "chrnos.h" 1309C 1310 LABINT = 'FC '//NAMN(ATMIND)(1:2) 1311 & //CHRNOS(ATMIND/100) 1312 & //CHRNOS(MOD(ATMIND,100)/10) 1313 & //CHRNOS(MOD(ATMIND,10)) 1314C 1315 RETURN 1316 END 1317C /* Deck SETDEG */ 1318 SUBROUTINE SETDEG 1319#include "implicit.h" 1320#include "priunit.h" 1321#include "maxaqn.h" 1322#include "maxmom.h" 1323#include "mxcent.h" 1324#include "maxorb.h" 1325#include "parnmr.h" 1326#include "nuclei.h" 1327#include "symmet.h" 1328#include "chrnos.h" 1329#include "chrxyz.h" 1330 1331C 1332 DO 50 I = 1, ESRNUC 1333 DEGEN(I) = 0.0D0 1334 DO 100 IREP=0,MAXREP 1335 ISCOR1 = IPTCNT(3*(NUCINF(I) - 1) + 1,IREP,2) 1336 IF (ISCOR1 .GT. 0) DEGEN(I)=DEGEN(I)+1.0D0 1337 100 CONTINUE 1338 50 CONTINUE 1339C 1340 RETURN 1341 END 1342C /* Deck SETISO */ 1343 SUBROUTINE SETISO 1344#include "implicit.h" 1345#include "priunit.h" 1346#include "parnmr.h" 1347#include "mxcent.h" 1348#include "nuclei.h" 1349C 1350 DO 100 I=1,ESRNUC 1351 IF (ISODAT(1,NUCINF(I)) .EQ. 0) THEN 1352 ISODAT(1,NUCINF(I))=1 1353 IATIS=1 1354 200 XATGV=DISOTP(INT(CHARGE(NUCINF(I))),IATIS,'GVAL') 1355 IF (DABS(XATGV) .LT. 1.0D-5) THEN 1356 IATIS=IATIS+1 1357 GO TO 200 1358 ENDIF 1359 ISODAT(2,NUCINF(I))=IATIS 1360 ENDIF 1361 100 CONTINUE 1362C 1363 RETURN 1364 END 1365C /* Deck SDOPER */ 1366 SUBROUTINE SDOPER(ATMIND,LABINT,NCOOR) 1367#include "implicit.h" 1368#include "priunit.h" 1369#include "maxaqn.h" 1370#include "maxmom.h" 1371#include "mxcent.h" 1372#include "maxorb.h" 1373#include "parnmr.h" 1374 INTEGER ATMIND, NCOOR 1375 CHARACTER LABINT(9*ATMNUM)*8 1376#include "nuclei.h" 1377#include "symmet.h" 1378#include "chrnos.h" 1379#include "chrxyz.h" 1380 1381C 1382 NCOOR=0 1383 DEGEN(ATMIND) = 0.0D0 1384 DO 100 IREP=0,MAXREP 1385 DO 200 ICOOR1=1,3 1386 ISCOR1 = IPTCNT(3*(NUCINF(ATMIND ) - 1) + ICOOR1,IREP,2) 1387 IF (ISCOR1 .GT. 0) THEN 1388 IF (ICOOR1 .EQ. 1) DEGEN(ATMIND)=DEGEN(ATMIND)+1.0D0 1389 DO 300 ICOOR2 = 1, 3 1390 NCOOR=NCOOR+1 1391 LABINT(NCOOR) = 'SD '//CHRNOS(ISCOR1/100) 1392 & //CHRNOS(MOD(ISCOR1,100)/10) 1393 & //CHRNOS(MOD(ISCOR1,10)) 1394 & //' '//CHRXYZ(-ICOOR2) 1395 300 CONTINUE 1396 SDIND(1,ISCOR1)=ATMIND 1397 SDIND(2,ISCOR1)=ICOOR1 1398 ENDIF 1399 200 CONTINUE 1400 100 CONTINUE 1401C 1402 RETURN 1403 END 1404 1405 1406C /* Deck pcmgdt */ 1407 SUBROUTINE PCMGDT(CREF,CMO,INDXCI,DV,G,ESOLT,WRK,LFREE, 1408 $ NHCREF,KREFSY) 1409C 1410C Based on SOLGDT: 1411C 13-02-2006 Luca Frediani 1412C 1413C Purpose: calculate MCSCF energy and gradient contribution 1414C from a surrounding medium with PCM 1415C 1416C Output: 1417C G MCSCF gradient with solvation contribution added 1418C 1419#include "implicit.h" 1420#include "dummy.h" 1421C 1422 DIMENSION CREF(*), CMO(*), INDXCI(*) 1423 DIMENSION DV(*), G(*), WRK(*) 1424 PARAMETER ( D1=1.0d0, DM1=-1.0d0, D0 = 0.0D0, D2 = 2.0D0 ) 1425#include "thrzer.h" 1426#include "iratdef.h" 1427#include "priunit.h" 1428#include "infrsp.h" 1429#include "mxcent.h" 1430#include "orgcom.h" 1431#include "pcmdef.h" 1432#include "pcmlog.h" 1433#include "pcm.h" 1434C 1435C Used from common blocks: 1436C INFINP: NLMSOL, LSOLMX, EPSOL, RSOL(3) 1437C INFVAR: NCONF, NWOPT, NVAR, NVARH 1438C INFORB: NNASHX, NNBASX, NNORBX, etc. 1439C INFIND: IROW(*) 1440C INFTAP: LUSOL, LUIT2 1441C INFPRI: IPRSOL 1442C 1443#include "maxash.h" 1444#include "maxorb.h" 1445#include "infinp.h" 1446#include "infvar.h" 1447#include "inforb.h" 1448#include "infind.h" 1449#include "inftap.h" 1450#include "infpri.h" 1451C 1452 LOGICAL FNDLAB, EXP1VL, TOFILE, TRIMAT 1453 PARAMETER (MXLMAX = 50) 1454 CHARACTER*8 STAR8, SOLVDI, EODATA, LABINT(9*MXCENT) 1455 DATA STAR8/'********'/ 1456 DATA SOLVDI/'SOLVDIAG'/, EODATA/'EODATA '/ 1457C 1458C Statement functions; 1459C define automatic arrays (dynamic core allocation) 1460C 1461 CALL QENTER('PCMGDT') 1462C 1463C 1464C Core allocation 1465C DIASH diagonal of solvent contribution to Hessian 1466C GRDLM TELM gradient for current l,m value in the l,m loop 1467C UCMO CMO unpacked (i.e. no symmetry blocking) 1468C 1469 KJENAO = 1 1470 KJEN = KJENAO + NNBASX 1471 KJ1AO = KJEN + NNORBX 1472 KJ1 = KJ1AO + NNBASX*NSYM 1473 KJENAC = KJ1 + NNORBX 1474 KJ1AC = KJENAC + NNASHX 1475 KDIASH = KJ1AC + NNASHX 1476 KGRDLM = KDIASH + NVAR 1477 KUCMO = KGRDLM + NVARH 1478 KJ2GRD = KUCMO + NORBT*NBAST 1479 KPOT = KJ2GRD + NVAR 1480 KDENC = KPOT + NTS 1481 KDENV = KDENC + N2BASX 1482 KW10 = KDENV + N2BASX 1483C 1.3 rest of CSF contribution 1484 KDIALM = KW10 1485 KTDV = KDIALM + NVAR 1486 KW20 = KTDV + NASHT * NASHT 1487C 1488C Allocations for non-equlibrium energy solvation 1489 IF (NONEQ) THEN 1490 KMPOT = KW20 1491 KQSEGR = KMPOT + NTS * NTS 1492 KPOTGR = KQSEGR + NTS 1493 KW21 = KPOTGR + NTS 1494 ELSE 1495 KW21 = KW20 1496 END IF 1497C 1498 LW21 = LFREE - KW21 1499C 1500 KWRK1 = KW21 1501 LWRK1 = LFREE - KWRK1 1502 IF (LWRK1 .LT. 0) CALL ERRWRK('PCMGDT',-KWRK1,LFREE) 1503C 1504 1505 KFREE=1 1506 ISPIN1=0 1507 ISPIN2=0 1508 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, 1509 * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 1510 * INDXCI,WRK(KWRK1),KFREE,LWRK1) 1511C 1512C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX 1513C 1514 CALL DSITSP(NASHT,WRK(KTDV),DV) 1515C 1516clf IF (IPRPCM .GE. 130) THEN 1517 IF (.true.) THEN 1518 WRITE (LUPRI,'(/A/A,2I10)') 1519 * ' --- PCMGDT - gtot (input) - non-zero elements', 1520 * ' NCONF, NWOPT =',NCONF,NWOPT 1521 DO 40 I = 1,NCONF 1522 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1523 * ' conf #',I,G(I) 1524 40 CONTINUE 1525 DO 50 I = NCONF+1,NVAR 1526 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1527 * ' orb #',I,G(I) 1528 50 CONTINUE 1529 END IF 1530 1531 IF (IPRSOL .GE. 20 .AND. NASHT .GT. 0) THEN 1532 WRITE (LUPRI,'(/A)') ' SOLGDT - DV matrix :' 1533 CALL OUTPAK(DV,NASHT,1,LUPRI) 1534 END IF 1535C 1536C Unpack symmetry blocked CMO 1537C Loop over l,m expansion 1538C 1539 CALL UPKCMO(CMO,WRK(KUCMO)) 1540 CALL DZERO(WRK(KDIASH),NVAR) 1541 CALL DZERO(WRK(KDIALM),NVAR) 1542 CALL DZERO(WRK(KJEN),NNORBX) 1543C 1544C Read JEN = J(en) + J(ne) in ao basis and transform to mo basis. 1545C Extract active-active block in WRK(KJENAC). 1546C 1547 LUPBKP = LUPROP 1548 IF (LUPROP .LT. 0) CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN', 1549 & 'SEQUENTIAL','UNFORMATTED',IDUMMY,.FALSE.) 1550CLF IF (.NOT. (FNDLAB('NE-PCMIN',LUPROP))) THEN 1551 IF (.TRUE.) THEN 1552 CALL PCMJMAT(WRK(KJENAO),NNBASX,WRK(KWRK1),LWRK1) 1553 END IF 1554 REWIND (LUPROP) 1555 CALL REAPCM('NE-PCMIN','PCMGRD ',LUPROP,WRK(KJENAO),NNBASX) 1556 CALL UTHU(WRK(KJENAO),WRK(KJEN),WRK(KUCMO),WRK(KWRK1), 1557 & NBAST,NORBT) 1558 IF (NASHT .GT. 0) THEN 1559 CALL GETAC2(WRK(KJEN),WRK(KJENAC)) 1560 END IF 1561 IF (IPRPCM .GE. 15) THEN 1562 WRITE (LUPRI,'(/A)') ' JEN_ao matrix:' 1563 CALL OUTPAK(WRK(KJENAO),NBAST,1,LUPRI) 1564 WRITE (LUPRI,'(/A)') ' JEN_mo matrix:' 1565 CALL OUTPAK(WRK(KJEN), NORBT,1,LUPRI) 1566 IF (NASHT .GT. 0) THEN 1567 WRITE (LUPRI,'(/A)') ' JEN_ac matrix:' 1568 CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) 1569 END IF 1570 END IF 1571C 1572C Expextation value of JEN (=PB) 1573C 1574 TJEN = SOLELM(DV,WRK(KJENAC),WRK(KJEN),TJENAC) 1575 IF (IPRPCM .GE. 6) THEN 1576 WRITE (LUPRI,'(A,F17.8)') 1577 * ' --- JEN expectation value(=PB) :',TJEN 1578 WRITE (LUPRI,'(A,F17.8)') 1579 * ' --- active part of JEN(=PB) :',TJENAC 1580 END IF 1581Cbm PB=-TJEN 1582Cbm WRITE(LUPRI,*)'PB =',PB 1583C 1584C Read J2 in ao basis and transform to mo basis. 1585C Extract active-active block in WRK(KJ2AC). 1586C 1587 CALL DZERO(WRK(KDENC),N2BASX) 1588 CALL DZERO(WRK(KDENV),N2BASX) 1589 CALL FCKDEN((NISHT.GT.0),(NASHT.GT.0),WRK(KDENC),WRK(KDENV), 1590 & CMO,DV,WRK(KWRK1),LW21) 1591 CALL DAXPY(N2BASX,1.0D0,WRK(KDENV),1,WRK(KDENC),1) 1592 CALL DZERO(WRK(KDENV),N2BASX) 1593 CALL DGEFSP(NBAST,WRK(KDENC),WRK(KDENV)) 1594ckr CALL DZERO(WRK(KDENC),N2BASX) 1595ckr CALL PKSYM1(WRK(KDENV),WRK(KDENC),NBAS,NSYM,1) 1596 1597 EXP1VL = .TRUE. 1598 TOFILE = .FALSE. 1599 CALL J1INT(WRK(KPOT),EXP1VL,WRK(KDENV),1,TOFILE,'NPETES ', 1600 & 1,WRK(KWRK1),LW21) 1601 CALL DZERO(QSE,MXTS) 1602 CALL V2Q(WRK(KWRK1),WRK(KPOT),QSE,QET,.FALSE.) 1603 CALL GPCLOSE(LUPCMD,'KEEP') 1604C CALL DSCAL(NTS,-1.0D0,QSE,1) 1605C 1606C Read J1 (=ELECTROSTATIC POTENTIAL) in ao basis and transform to mo 1607C basis. Extract active-active block in WRK(KJ1AC). 1608C 1609 XI = DIPORG(1) 1610 YI = DIPORG(2) 1611 ZI = DIPORG(3) 1612Ckr 1613Ckr should be parallelized, but a lot of information to send 1614Ckr However, in general not used for SCF and DFT optimizations 1615Ckr 1616 DO I = 1 , NTSIRR 1617 L = 1 1618 NCOMP = NSYM 1619 DIPORG(1) = XTSCOR(I) 1620 DIPORG(2) = YTSCOR(I) 1621 DIPORG(3) = ZTSCOR(I) 1622 EXP1VL = .FALSE. 1623 TOFILE = .FALSE. 1624 KPATOM = 0 1625 TRIMAT = .TRUE. 1626 CALL GET1IN(WRK(KJ1AO),'NPETES ',NCOMP,WRK(KWRK1),LW21,LABINT, 1627 & INTREP,INTADR,L,TOFILE,KPATOM,TRIMAT,DUMMY,EXP1VL, 1628 & DUMMY,IPRPCM) 1629 CALL UTHU(WRK(KJ1AO),WRK(KJ1),WRK(KUCMO),WRK(KWRK1), 1630 & NBAST,NORBT) 1631 IF (NASHT .GT. 0) THEN 1632 CALL GETAC2(WRK(KJ1),WRK(KJ1AC)) 1633 END IF 1634 IF (NONEQ) THEN 1635Clf here we will need the singlet density, I beleive..... 1636 WRK(KPOT + I - 1) = SOLELM(DV,WRK(KJ1AC),WRK(KJ1),TJ1AC) 1637 END IF 1638 1639 IF (IPRPCM .GE. 15) THEN 1640 WRITE (LUPRI,'(/A)') ' J1_ao matrix:' 1641 CALL OUTPAK(WRK(KJ1AO),NBAST,1,LUPRI) 1642 WRITE (LUPRI,'(/A)') ' J1_mo matrix:' 1643 CALL OUTPAK(WRK(KJ1), NORBT,1,LUPRI) 1644 IF (NASHT .GT. 0) THEN 1645 WRITE (LUPRI,'(/A)') ' J1_ac matrix:' 1646 CALL OUTPAK(WRK(KJ1AC),NASHT,1,LUPRI) 1647 END IF 1648 END IF 1649 CALL DAXPY(NNORBX,-QSE(I),WRK(KJ1),1,WRK(KJEN),1) 1650C 1651C Due to the computational cost of building the CI gradient, and the lack 1652C of importance on convergence rates, we skip the construction of the 1653C solvent contribution to the diagonal Hessian. K.Ruud, Oct.-01 1654C 1655C CALL DZERO(WRK(KGRDLM),NVARH) 1656C IF (NCONF .GT. 1) THEN 1657C CALL SOLGC(CREF,WRK(KJ1AC),TJ1AC,WRK(KGRDLM),INDXCI, 1658C & WRK(KWRK1),LWRK1) 1659CC CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) 1660CC END IF 1661C IF (NWOPT .GT. 0) THEN 1662C CALL SOLGO(DCVAL,DV,WRK(KJ1),WRK(KGRDLM+NCONF)) 1663C END IF 1664C READ(LUGRDQ)WRK(KJ2GRD) 1665C DO J = NCONF, NVAR - 1 1666C WRK(KDIASH+J) = WRK(KDIASH+J) - WRK(KGRDLM+J)*WRK(KJ2GRD+J) 1667C ENDDO 1668 ENDDO 1669Ckr 1670Ckr End of parallelization loop 1671Ckr 1672C 1673C Expextation value of J + X(0) = PB + PX 1674C 1675 IF (NASHT .GT. 0) THEN 1676 CALL GETAC2(WRK(KJEN),WRK(KJENAC)) 1677 END IF 1678 TJEN = SOLELM(DV,WRK(KJENAC),WRK(KJEN),TJENAC) 1679 1680 IF (IPRPCM .GE. 6) THEN 1681 CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) 1682 WRITE (LUPRI,'(A,F17.8)') 1683 * ' --- TJEN expectation value :',TJEN 1684 WRITE (LUPRI,'(A,F17.8)') 1685 * ' --- active part of TJEN :',TJENAC 1686 END IF 1687C 1688 KFREE=1 1689 ISPIN1=1 1690 ISPIN2=0 1691 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, 1692 * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., 1693 * INDXCI,WRK(KWRK1),KFREE,LWRK1) 1694Clf 1695 IF (.true.) THEN 1696 WRITE (LUPRI,'(/A)') ' JEN_ao matrix:' 1697 CALL OUTPAK(WRK(KJENAO),NBAST,1,LUPRI) 1698 WRITE (LUPRI,'(/A)') ' JEN_mo matrix:' 1699 CALL OUTPAK(WRK(KJEN), NORBT,1,LUPRI) 1700 IF (NASHT .GT. 0) THEN 1701 WRITE (LUPRI,'(/A)') ' JEN_ac matrix:' 1702 CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) 1703 END IF 1704 END IF 1705 1706 1707C 1708C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX 1709C 1710 CALL DSITSP(NASHT,WRK(KTDV),DV) 1711 CALL DZERO(WRK(KGRDLM),NVARH) 1712 IF (NCONF .GT. 1) THEN 1713 CALL SOLGC(CREF,WRK(KJENAC),TJENAC,WRK(KGRDLM),INDXCI, 1714 & WRK(KWRK1),LWRK1) 1715C CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) 1716 END IF 1717 IF (NWOPT .GT. 0) THEN 1718 CALL SOLGO(D0,DV,WRK(KJEN),WRK(KGRDLM+NCONF)) 1719 END IF 1720C 1721C 1722C Obtain DIALM = diagonal TE(l,m) Hessian 1723C = 2 ( <i|R(l,m)|i> - TE(l,m) ) 1724C Add the DIALM contribution and the GRDLM contribution 1725C to solvent Hessian diagonal. 1726C 1727C CALL SOLDIA(TJENAC,WRK(KJENAC),INDXCI, 1728C * WRK(KJEN),DV,WRK(KDIALM),WRK(KW21),LW21) 1729C CALL SOLDIA(TELM,RLMAC,INDXCI,RLM,DV,DIAG,WRK,LFREE) 1730C 1731 DO 420 I = 0,(NVAR-1) 1732 WRK(KDIASH+I) = WRK(KDIASH+I) 1733 * - WRK(KDIALM+I) 1734 420 CONTINUE 1735C 1736C test orthogonality 1737C 1738 IF (IPRPCM .GE. 120) THEN 1739 WRITE (LUPRI,'(/A)')' --- PCMGRD - grdj1, grdj2, dialm, '// 1740 & 'diash, cref' 1741 DO 430 I = 1,NCONF 1742 WRITE (LUPRI,'(A,I10,3F10.6)') ' conf #',I, 1743 * WRK(KDIALM-1+I), 1744 * WRK(KDIASH-1+I),CREF(I) 1745 430 CONTINUE 1746 END IF 1747C 1748 TEST = DDOT(NCONF,CREF,1,WRK(KGRDLM),1) 1749 IF (ABS(TEST) .GT. 1.D-8) THEN 1750 NWARN = NWARN + 1 1751 WRITE (LUPRI,'(/A,/A,1P,D12.4)') 1752 * ' --- PCMGRD WARNING, for B', 1753 * ' <CREF | GRADB > =',TEST 1754 END IF 1755C 1756C TEST = DDOT(NCONF,CREF,1,WRK(KGRDJ1),1) 1757 IF (ABS(TEST) .GT. 1.D-8) THEN 1758 NWARN = NWARN + 1 1759 WRITE (LUPRI,'(/A,/A,1P,D12.4)') 1760 * ' --- PCMGRD WARNING, for J1 ', 1761 * ' <CREF | GRADJ1 > =',TEST 1762 END IF 1763C TEST = DDOT(NCONF,CREF,1,WRK(KGRDJ2),1) 1764 IF (ABS(TEST) .GT. 1.D-8) THEN 1765 NWARN = NWARN + 1 1766 WRITE (LUPRI,'(/A,/A,1P,D12.4)') 1767 * ' --- PCMGRD WARNING, for J2', 1768 * ' <CREF | GRADJ2 > =',TEST 1769 END IF 1770C 1771C Add PCM gradient contribution to MCSCF gradient 1772C 1773 CALL DAXPY(NVARH,DM1,WRK(KGRDLM),1,G,1) 1774clf IF (IPRPCM .GE. 140) THEN 1775 IF (.true.) THEN 1776 WRITE (LUPRI,'(/A/A,2I10)') 1777 * ' --- PCMGRD - grdB, gtot (accum) - non-zero grdlm', 1778 * ' NCONF, NWOPT =',NCONF,NWOPT 1779 DO 440 I = 1,NCONF 1780 IF (WRK(KGRDLM-1+I) .NE. D0) 1781 * WRITE (LUPRI,'(A,I10,3F15.10)') 1782 * ' conf #',I,WRK(KGRDLM-1+I),G(I) 1783 440 CONTINUE 1784 DO 450 I = NCONF+1,NVAR 1785 IF (WRK(KGRDLM-1+I) .NE. D0) 1786 * WRITE (LUPRI,'(A,I10,3F15.10)') 1787 * ' orb #',I,WRK(KGRDLM-1+I),G(I) 1788 450 CONTINUE 1789 END IF 1790C 1791C 1792C 1793 IF (IPRPCM .GE. 130) THEN 1794 WRITE (LUPRI,'(/A/A,2I10)') 1795 * ' --- PCMGRD - gtot (output) - non-zero elements', 1796 * ' NCONF, NWOPT =',NCONF,NWOPT 1797 DO 840 I = 1,NCONF 1798 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1799 * ' conf #',I,G(I) 1800 840 CONTINUE 1801 DO 850 I = NCONF+1,NVAR 1802 IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') 1803 * ' orb #',I,G(I) 1804 850 CONTINUE 1805 END IF 1806 IF (LUPBKP .LT. 0) CALL GPCLOSE(LUPROP,'KEEP') 1807 CALL QEXIT('PCMGDT') 1808 RETURN 1809C end of pcmgdt. 1810 END 1811! --- end of DALTON/rsp/lagrang.F --- 1812