! ! Dalton, a molecular electronic structure program ! Copyright (C) by the authors of Dalton. ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU Lesser General Public ! License version 2.1 as published by the Free Software Foundation. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! Lesser General Public License for more details. ! ! If a copy of the GNU LGPL v2.1 was not distributed with this ! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. ! ! ! ! File: dalton/rsp/lagrang.F ! ! The routines in this file are generally related to ! setup the Lagrangian for doublet reference states with ! triplet property operators, and solving for the corresponding ! response vector. ! ! This is what has been called the "restricted-unrestricted" approach ! iin the literature. ! ! The driver RSPESR uses this for calculating expectation vaules of triplet ! operators for non-singlet reference states. ! This file also includes the input routine ESRINP for specification ! of ESR calculations, using RSPESR here for the requested expectation values, ! as well as routines in rspg.F and in rspzfs.F for g-tensor and ZFS, ! depending on input. ! ! C=========================================================================== CRevision 1.4 2000/05/24 18:45:09 hjj Cnew getref calls with appropriate NDREF (fixing CSF and triplet) Cstop if solvent and CSF Cbugfix: too little was allocated for KGP2 C=========================================================================== C C C /* Deck lagran */ SUBROUTINE LAGRAN(WORD,FC,FV,CMO,UDV,PV,XINDX,WRK,LWRK) C C 18-AUG-1991 C C Called from GETGPV when label(5:8) = 'LAGR' C #include "implicit.h" #include "dummy.h" DIMENSION CMO(*),UDV(NASHT,*),FC(*),FV(*),PV(*),XINDX(*),WRK(*) CHARACTER*8 WORD C-- common blocks: #include "maxorb.h" #include "priunit.h" #include "infdim.h" #include "infvar.h" #include "inforb.h" #include "infrsp.h" #include "wrkrsp.h" #include "infpri.h" #include "infinp.h" #include "dftcom.h" #include "pcmlog.h" C PARAMETER (DHALF=0.5D0,D1=1.0D0,DM1=-1.0D0) C CALL QENTER('LAGRAN') C IF (IPRRSP .GE. 10) WRITE (LUPRI,7010) WORD 7010 FORMAT(//' Output from LAGRAN:'/' ===================' * //' Property label = ',A/) IF (SOPPA) THEN WRITE (LUPRI,*) 'SOPPA not implemented in LAGRAN yet, sorry!' CALL QUIT('LAGRAN-ERROR: SOPPA not implemented') END IF clf IF ((FLAG(16).OR.PCM).AND.(NASHT.GT.1).AND.(NCREF.NE.KZCONF)) THEN Chj: see comments in RSPSLV in rspsol.F clf WRITE(LUPRI,*)'Solvent not implemented in LAGRAN with CSFs' clf WRITE(LUPRI,*)'Use .DETERMINANTS and try again!' clf CALL QUIT('Solvent not implemented in LAGRAN with CSFs') clf END IF C C ALLOCATE WORK SPACE C KGP = 1 KCREF = KGP + KZYVAR NDREF = MAX(KZCONF,NCREF) Chj ... RSPOLI requires determinants for ZYCVEC when triplet, Chj ... thus KZCONF. However, for ROHF NCREF = 1, KZCONF = 0 Chj and we need CREF(1) for RSPDM. Thus MAX(KZCONF,NCREF) KTOT = KCREF + NDREF C CALL DZERO(WRK(KGP),KZYVAR) CALL GETREF(WRK(KCREF),NDREF) C C CONSTRUCT the Lagrangian vector, the GP vector, C using routine for the ORBITAL PART OF LINEAR TRANSFORMED VECTORS C IF (.NOT.RSPCI) THEN C C ALLOCATE WORK SPACE FOR RSPOLI C NCSIM = 1, NOSIM = 0 C KFVTD = KTOT LFVTD = N2ORBX IF (DOMCSRDFT) LFVTD = 2*LFVTD C ... Extra allocation for "FCTD" in MCSCF-SRDFT, C is needed for the VxcTD matrix. KQATD = KFVTD + LFVTD KQBTD = KQATD + NORBT * NASHT KWRK1 = KQBTD + NORBT * NASHT LWRK1 = LWRK - KWRK1 IF (LWRK1.LT.0) CALL ERRWRK('LAGRAN',KWRK1-1,LWRK) IF (DOHFSRDFT .OR. DOMCSRDFT) THEN CALL QUIT('LAGRANGE is not implemented yet for srDFT') END IF IF ((DFTADD.OR.HSROHF).AND.TRPLET) THEN CALL DAXPY(NNORBT,D1,FC,1,FV,1) CALL DSCAL(NNORBT,DM1,FV,1) END IF CALL RSPOLI(1,0,UDV,WRK(KCREF),NDREF,FC,FV,PV,DUMMY, * DUMMY,DUMMY, DUMMY,DUMMY,WRK(KFVTD), * WRK(KQATD),WRK(KQBTD),WRK(KGP), * XINDX,CMO,WRK(KWRK1),LWRK1) C CALL RSPOLI(NCSIM,NOSIM,UDV,ZYCVEC,LZYCVEC,FC,FV,PVX,ZYMAT, C * FCONE,FVONE,QAONE,QBONE,FVTD,QATD,QBTD,EVECS, C * XINDX,CMO,WRK,LWRK) IF ((DFTADD.OR.HSROHF).AND.TRPLET) THEN CALL DAXPY(NNORBT,D1,FC,1,FV,1) CALL DSCAL(NNORBT,DM1,FV,1) END IF C END IF IF (IPRRSP.GT.120) THEN WRITE(LUPRI,'(/A)') ' LAGRAN: LAGRANGIAN VECTOR ' CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF C C section for calculating esr-properties with solvent C contributions. C IF (FLAG(16).OR.PCM) THEN C C workspace allocation C KGP2 = KTOT KTDV = KGP2 + NVARH Chj KTDV = KGP2 + KZVAR Chj ... SOLGDT uses NVARH (.gt. KZVAR !) for GP2 KDV = KTDV + N2ASHX KWRK1 = KDV + NNASHX IF (FLAG(16)) THEN KWRK2 = KWRK1 + NLMSOL * 2 ELSE IF (PCM) THEN KWRK2 = KWRK1 ELSE KWRK2 = KWRK1 END IF LWRK1 = LWRK - KWRK2 C C initialize solvent gradient. C CALL DZERO(WRK(KGP2),NVARH) C C ISPIN1 = 1 ISPIN2 = 0 CALL RSPDM(IREFSY,IREFSY,NDREF,NDREF,WRK(KCREF),WRK(KCREF), * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * XINDX,WRK,KWRK2,LWRK1) C CALL RSPDM(ISYM,ISYM,NCDIM.NCDIM,REF,REF, C TUDV,DUMMY, ISPIN1,ISPIN2,TDM,NORHO2, C XINDX,WRK,KFREE,LFREE) IF (IPRRSP.GT.110) THEN WRITE(LUPRI,'(/A)')'LAGRANG: ONE ELECTRON SPIN DENSITY' CALL OUTPUT(WRK(KTDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) END IF C C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX C CALL DGETSP(NASHT,WRK(KTDV),WRK(KDV)) C IF (IPRRSP.GT.50 .AND. NASHT.GT.0) THEN WRITE(LUPRI,'(/A)') & 'LAGRANG solvent: PACKED ONE ELECTRON SPIN DENSITY MATRIX' CALL OUTPAK(WRK(KDV),NASHT,1,LUPRI) END IF C C use wavefunction routine to calculate solvent gradient. C IF (FLAG(16)) THEN CALL SOLGDT(WRK(KCREF),CMO,XINDX,WRK(KDV),WRK(KGP2), * ESOLT,WRK(KWRK1),WRK(KWRK2),LWRK1,NDREF,IREFSY) ELSE IF (PCM) THEN CALL PCMGDT(WRK(KCREF),CMO,XINDX,WRK(KDV),WRK(KGP2), * ESOLT,WRK(KWRK2),LWRK1,NDREF,IREFSY) ELSE CALL QUIT("Lagran: Either FLAG(16) or PCM must be TRUE") END IF C C scale by a half for consistency with response and add to GP. C CALL DAXPY(KZWOPT,DHALF,WRK(KGP2+KZCONF),1,WRK(KGP+KZCONF),1) CALL DAXPY(KZWOPT,(-DHALF),WRK(KGP2+KZCONF),1, * WRK(KGP+KZVAR+KZCONF),1) C IF (IPRRSP.GT.120) THEN WRITE(LUPRI,'(/A)') ' LAGRAN: 2 * SOLVENT LAGRANGIAN VECTOR' CALL OUTPUT(WRK(KGP2),1,NVARH,1,1,NVARH,1,1,LUPRI) WRITE(LUPRI,'(/A)') & ' LAGRAN: Electronic + solvent LAGRANGIAN VECTOR' CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF C ENDIF C C *** END OF LAGRAN C CALL QEXIT('LAGRAN') RETURN END C /* Deck rsplan */ SUBROUTINE RSPLAN(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) C C Lagrangian solution vector is returned in WRK(1) C #include "implicit.h" #include "dummy.h" #include "iratdef.h" #include "inftap.h" C DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) DIMENSION XINDX(*),WRK(*) C #include "priunit.h" #include "pgroup.h" #include "infpri.h" #include "infrsp.h" #include "wrkrsp.h" #include "rspprp.h" #include "inflr.h" #include "infdim.h" #include "inforb.h" #include "infesr.h" C C Local variables C CHARACTER*8 LABEL, BLANK PARAMETER ( D0 = 0.0D0 , DM8 = 1.0D-8, BLANK = ' ' ) C CALL QENTER('RSPLAN') C C DETERMINE SECOND ORDER MOLECULAR PROPERTIES C KREDE = 1 KREDS = KREDE + MAXRM*MAXRM KIBTYP = KREDS + MAXRM*MAXRM KEIVAL = KIBTYP + MAXRM KRESID = KEIVAL + MAXRM KEIVEC = KRESID + MAXRM KREDGP = KEIVEC + MAXRM*MAXRM KGP = KREDGP + MAXRM KWRK1 = KGP + KZYVAR LWRK1 = LWRK - KWRK1 C IF (LWRK1.LT.3*KZYVAR) THEN WRITE (LUERR,9100) LWRK1,3*KZYVAR CALL QTRACE(LUERR) CALL QUIT('RSPLAN ERROR, INSUFFICIENT MEMORY') ENDIF C C WORK SPACE FOR RSPEVE C KWRKE = KWRK1 KBVECS = KWRKE + KZYVAR 9100 FORMAT(/' RSPLAN, work space too small for 3 (Z,Y)-vectors', * /' had:',I10,', need more than:',I10) C KZRED = 0 KZYRED = 0 THCRSP = THCESR IPRRSP = IPRESR MAXIT = MAXESR C C Call RSPCTL to solve linear set of response equations C IF (TRPLET) THEN LABEL = 'TRIPLAGR' ELSE LABEL = 'SINGLAGR' END IF WRITE (LUPRI,'(//A,I2,3A/2A)') & ' RSPLAN -- linear response calculation for symmetry',KSYMOP, & ' ( ',REP(KSYMOP-1),')', & ' RSPLAN -- operator label : ',LABEL CALL GETGPV(LABEL,FC,FV,CMO,UDV,PV,XINDX,ANTSYM,WRK(KGP),LWRK1) IF (IPRRSP.GT.120) THEN WRITE(LUPRI,'(/A)') ' RSPLAN: LAGRANGIAN VECTOR ' CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF XNOR = DDOT(KZYVAR,WRK(KGP),1,WRK(KGP),1) IF ( XNOR.LT.DM8 ) THEN WRITE(LUPRI,'(/A)') * ' LAGRANGIAN VECTOR IS ZERO, SOLUTION IS SET TO ZERO' CALL DZERO(WRK(1),KZYVAR) GO TO 9999 END IF WRK(KEIVAL) = D0 KEXSIM = 1 KEXCNV = 1 CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, * .TRUE.,LABEL,BLANK,WRK(KGP),WRK(KREDGP), * WRK(KREDE),WRK(KREDS), * WRK(KIBTYP),WRK(KEIVAL),WRK(KRESID),WRK(KEIVEC), * XINDX,WRK(KWRK1),LWRK1) C CALL RSPCTL(CMO,UDV,PV,FC,FV,FCAC,H2AC, C * LINEQ,GP,REDGP,REDE,REDS, C * IBTYP,EIVAL,EIVEC,XINDX,WRK,LWRK) C NBX = 1 IBOFF = 0 CALL RSPEVE(WRK(KIBTYP),WRK(KEIVAL),WRK(KEIVEC), * WRK(KBVECS),WRK(KWRKE),NBX,IBOFF) C CALL RSPEVE(IBTYP,EIVAL,EIVEC,BVECS,WRK,NBX,IBOFF) IF ( IPRRSP.GE.10 ) THEN WRITE (LUPRI,'(/A)') ' SOLUTION VECTOR FOR LAGRANGIAN ' CALL RSPPRC(WRK(KBVECS),KZCONF,KZVAR,LUPRI) CALL RSPPRO(WRK(KBVECS+KZCONF),KZVAR,UDV,LUPRI) END IF CALL DCOPY(KZYVAR,WRK(KBVECS),1,WRK(1),1) IF ( IPRRSP.GE.10 ) THEN CALL HEADER(LABEL, LUPRI) CALL OUTPUT(WRK, 1, KZVAR, 1, 2, KZVAR, 2, 1, LUPRI) END IF CALL WRTRSP( & & LURSP, KZYVAR, WRK, LABEL, ' ', D0, D0, 1, 1, D0, .TRUE.& &) C C *** end of RSPLAN -- C 9999 CALL QEXIT('RSPLAN') RETURN END C /* Deck esrinp */ SUBROUTINE ESRINP(WORD) C C Module for reading "*ESR " input under **RESPONSE C #include "implicit.h" C #include "priunit.h" #include "gnrinf.h" #include "infrsp.h" #include "rspprp.h" #include "inflr.h" #include "infesr.h" #include "infpri.h" #include "zfs.h" #include "gtensor.h" #include "maxorb.h" #include "infinp.h" #include "parnmr.h" C LOGICAL NEWDEF PARAMETER ( NTABLE = 10 ) CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 CHARACTER*8 LABEL CHARACTER LABDAT(9*ATMNUM)*8 C DATA TABLE /'.TRPPRP', '.SNGPRP', '.MAX IT', '.THCESR', * '.PRINT ', '.G-TENS', '.ZFS ', '.FCCALC', * '.ATOMS ','.SDCALC' / C C Initialize new HFC input FCFLG = .FALSE. SDFLG = .FALSE. ESRNUC = 0 DO 50 I = 1, ATMNUM DO 55 IH = 1, 7 ISODAT(IH,I) = 0 55 CONTINUE 50 CONTINUE C NEWDEF = (WORD .EQ. '*ESR ') ICHANG = 0 ESRCAL = .FALSE. GCALC = .FALSE. ZFSCAL = .FALSE. IF (NEWDEF) THEN WORD1 = WORD 100 CONTINUE READ (LUCMD, '(A7)') WORD CALL UPCASE(WORD) PROMPT = WORD(1:1) IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') GO TO 100 IF (PROMPT .EQ. '.') THEN ICHANG = ICHANG + 1 DO 200 I = 1, NTABLE IF (TABLE(I) .EQ. WORD) THEN GO TO (1,2,3,4,5,6,7,8,9,10), I END IF 200 CONTINUE IF (WORD .EQ. '.OPTION') THEN CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) GO TO 100 END IF WRITE (LUPRI,'(/,3A,/)') ' KEYWORD "',WORD, * '" NOT RECOGNIZED IN ESRINP.' CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) CALL QUIT(' ILLEGAL KEYWORD IN LRINP ') 1 CONTINUE ! .TRPPRP READ(LUCMD,'( A )')LABEL ESROPT( INDPRP(LABEL)) = .TRUE. GO TO 100 2 CONTINUE ! .SNGPRP READ(LUCMD,'( A )')LABEL ESROPS( INDPRP(LABEL)) = .TRUE. GO TO 100 3 CONTINUE ! .MAX IT READ (LUCMD,*) MAXESR GO TO 100 4 CONTINUE ! .THCESR READ (LUCMD,*) THCESR GO TO 100 5 CONTINUE ! .PRINT READ (LUCMD,*) IPRESR GO TO 100 6 CONTINUE ! .G-TENSOR CALL GINP(WORD) GO TO 100 7 CONTINUE ! .ZFS ZFSCAL = .TRUE. GO TO 100 8 CONTINUE ! .FCCALC FCFLG = .TRUE. GO TO 100 9 CONTINUE ! .ATOMS READ (LUCMD,*) ESRNUC IF (ESRNUC .GT. ATMNUM) THEN WRITE(LUPRI,'(/A,I4,A,I4/A)') & ' Sorry, but .ATOMS under *ESR is',ESRNUC, & ' which is bigger than ATMNUM in parnmr.h',ATMNUM, & ' Either reduce .ATOMS or increase ATMNUM and recompile.' CALL QUIT('.ATOMS too big, see output') END IF READ (LUCMD,*) (NUCINF(IG), IG = 1, ESRNUC) GO TO 100 10 CONTINUE ! .SDCALC SDFLG = .TRUE. GO TO 100 ELSE IF (PROMPT .EQ. '*') THEN GO TO 300 ELSE WRITE (LUPRI,'(/3A/)') ' PROMPT "',WORD, * '" NOT RECOGNIZED IN ESRINP FOR *ESR.' CALL QUIT(' ILLEGAL PROMPT under *ESR') END IF GO TO 100 END IF 300 CONTINUE IF (THR_REDFAC .GT. 0.0D0) THEN ICHANG = ICHANG + 1 WRITE (LUPRI,'(3A,1P,D10.2)') '@ INFO ',WORD1, & ' thresholds multiplied with general factor',THR_REDFAC THCESR = THCESR*THR_REDFAC END IF C C FC labels generation C IF (FCFLG) THEN DO I = 1, ESRNUC CALL FCOPER(NUCINF(I),LABEL) ESROPT( INDPRP(LABEL)) = .TRUE. END DO IF (.NOT. SDFLG) CALL SETDEG END IF C C SD labels generation C IF (SDFLG) THEN DO I = 1, ESRNUC CALL SDOPER(NUCINF(I),LABDAT,INUM) DO IG = 1, INUM ESROPT( INDPRP(LABDAT(IG))) = .TRUE. END DO END DO END IF C C Check degen array C DO I = 1, ESRNUC IF (DEGEN(I) .LT. 1.0) DEGEN(I)=1.0 END DO C C Setting isotopes data for HFC/pNMR printing C CALL SETISO C WRITE(LUPRI,'(/A)')' ********* ESRINP ********' C C NTOESR = 0 DO I = 1,NPRLBL IF (ESROPS(I)) NTOESR = NTOESR + 1 IF (ESROPT(I)) NTOESR = NTOESR + 1 END DO ESRCAL = NTOESR .GT. 0 .OR. FCFLG .OR. SDFLG IF (ESRCAL .OR. GCALC .OR. ZFSCAL) THEN CALL HEADER('Changes of defaults under *ESR :',0) IF ( ESRCAL ) THEN WRITE (LUPRI,'(/A)') * ' ESR - hyperfine calculation carried out' IF (MAXESR.NE.60) WRITE(LUPRI,'(/A,I6)') * ' MAXIMUM NUMBER OF ITERATIONS. MAXESR =',MAXESR IF (THCESR.NE.1.0D-5) WRITE(LUPRI,'(/A,D13.6)') * ' THRESHOLD FOR CONVERGENCE. THCESR =',THCESR IF (IPRESR.NE.2) WRITE(LUPRI,'(/A,I5)') * ' PRINT LEVEL IN ESR ROUTINES. IPRESR =',IPRESR END IF IF (ZFSCAL) THEN WRITE (LUPRI,'(/A)') * ' ESR - zero-field splitting calculation carried out' END IF IF (GCALC) THEN WRITE (LUPRI,'(/A)') * ' ESR - g-tensor calculation carried out' IF (ECC) WRITE(LUPRI,'(A)') * ' Electronic charge centroid used as gauge origin' IF (MNFSO) WRITE(LUPRI,'(A)') * ' Mean-field spin-orbit approximation is used' IF (SCALED_CHARGES) WRITE(LUPRI,'(A)') * ' Two-electron spin-orbit by scaled nuclear charges' IF (.NOT.DOALL) THEN WRITE(LUPRI,'(/A)') * ' Selected contributions to electronic g_tensor:' IF (DORMC) WRITE(LUPRI,'(/A)') * ' Relativistic mass correction' IF (DOGC1) WRITE(LUPRI,'(/A)') * ' One-electron gauge correction' IF (DOGC2) WRITE(LUPRI,'(/A)') * ' Two-electron gauge correction' IF (DOOZSO1) WRITE(LUPRI,'(/A)') * ' One-electron spin-orbit + orbital Zeeman correction' IF (DOOZSO2) WRITE(LUPRI,'(/A)') * ' Two-electron spin-orbit + orbital Zeeman correction' END IF END IF ELSE WRITE (LUPRI,'(//A/)') * 'INFO: "*ESR " input ignored because no operators requested.' END IF IF (ISPIN .EQ. 1) THEN NWARN = NWARN + 1 WRITE (LUPRI,'(//A/)') * ' WARNING: "*ESR " input ignored because'// & ' singlet reference state has no ESR properties !!!' ESRCAL = .FALSE. FCFLG = .FALSE. SDFLG = .FALSE. GCALC = .FALSE. ZFSCAL = .FALSE. END IF C C *** END OF ESRINP C RETURN END C /* Deck rspesr */ SUBROUTINE RSPESR(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) C C Calculate triplet expectation values with Lagrangian correction C and singlet expectation values for open shell systems with MCSCF or CI C (non-zero spin densities). C C Revised March 2003 hjaaj C #include "implicit.h" #include "iratdef.h" C DIMENSION CMO(*),UDV(*),PV(*),FC(*),FV(*),FCAC(*),H2AC(*) DIMENSION XINDX(*),WRK(*) C PARAMETER ( D0 = 0.0D0, DUMMY = 1.0D+20 ) LOGICAL TRPSAVE C #include "maxorb.h" #include "infinp.h" C #include "priunit.h" #include "infpri.h" #include "infrsp.h" #include "wrkrsp.h" #include "rspprp.h" #include "inflr.h" #include "infdim.h" #include "inforb.h" #include "infesr.h" C C C Common blocks for HFC/pNMR printing C #include "codata.h" #include "gfac.h" #include "mxcent.h" #include "nuclei.h" #include "chrxyz.h" #include "parnmr.h" C C Local variables for HFC values handling C CHARACTER*8 LABEL CHARACTER*2 CTMP C CALL QENTER('RSPESR') CALL HEADER('Output from RSPESR module',0) C HFC printing counters IFCIND = 0 REFSPIN=DBLE(ISPIN-1) C C Zero SDVAL C DO 50 I = 1, ATMNUM DO 55 IH = 1, 3 DO 60 IG = 1, 3 SDVAL(IG,IH,I) = 0.0D0 60 CONTINUE 55 CONTINUE 50 CONTINUE C C DETERMINE SECOND ORDER MOLECULAR PROPERTIES C IF ( NESRT(KSYMOP).GT. 0 ) THEN NDREF = KZCONF C ... we use determinants when triplet, thus not NCREF NDREF = MAX(KZCONF,NCREF) Chj ... we use determinants triplet, thus KZCONF, Chj ... However, for ROHF NCREF = 1, KZCONF = 0 and Chj we need CREF(1) for RSPDM. Thus MAX(KZCONF,NCREF) C C Note that GETGPV uses WRK(KGP) as scratch, thus KGP last allocation ! C IF (RSPCI) THEN KTUDV = 1 KCREF = KTUDV + N2ASHX KWRK1 = KCREF + NDREF KLAGR = KCREF KGP = KCREF ELSE KLAGR = 1 CALL RSPLAN(CMO,UDV,PV,FC,FV,FCAC,H2AC,XINDX,WRK,LWRK) IF (IPRRSP.GT.120) THEN WRITE(LUPRI,'(/A)') & 'RSPESR: solution vector for lagrangian :' CALL OUTPUT(WRK,1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF KTUDV = KLAGR + KZYVAR KCREF = KTUDV + N2ASHX KGP = KCREF KWRK1 = KGP + MAX(KZYVAR,NDREF) END IF LWRK1 = LWRK - KWRK1 IF (LWRK1.LT.0) CALL ERRWRK('RSPESR',KWRK1-1,LWRK) C C Get one electron spin density (triplet MS=0 density) C CALL GETREF(WRK(KCREF),NDREF) KFREE = 1 LFREE = LWRK1 ISPIN1 = 1 ISPIN2 = 0 CALL RSPDM(IREFSY,IREFSY,NDREF,NDREF,WRK(KCREF),WRK(KCREF), * WRK(KTUDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * XINDX,WRK(KWRK1),KFREE,LFREE) C CALL RSPDM(ISYM,ISYM,NCDIM.NCDIM,REF,REF, C TUDV,DUMMY, ISPIN1,ISPIN2,TDM,NORHO2, C XINDX,WRK,KFREE,LFREE) IF (IPRRSP.GE.5) THEN WRITE(LUPRI,'(/A)')'RSPESR: one electron spin density ' CALL OUTPUT(WRK(KTUDV),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) END IF TRPSAVE= TRPLET TRPLET = .TRUE. DO 100 IOP = 1,NESRT(KSYMOP) IF (.NOT. RSPCI ) THEN IF (IPRRSP.GT.150) THEN WRITE(LUPRI,'(/A)') * 'RSPESR: Lagrangian vector before product' CALL OUTPUT(WRK(KLAGR),1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF CALL GETGPV(LBESRT(KSYMOP,IOP),FC,FV,CMO,UDV,PV,XINDX, * ANTSYM,WRK(KGP),LWRK1 ) C CALL GETGPV(LABELOP,FC,CMO,UDV,PV,XINDX,ANTSYM,WRK,LWRK) IF (IPRRSP.GT.120) THEN WRITE(LUPRI,'(/2A)') * 'RSPESR: GP vector with label: ',LBESRT(KSYMOP,IOP) CALL OUTPUT(WRK(KGP),1,KZVAR,1,2,KZVAR,2,1,LUPRI) END IF FOPLG = -DDOT(KZYVAR,WRK(1),1,WRK(KGP),1) IF (IPRRSP.GT.5) WRITE(LUPRI,'(/A,1P,D18.10)') & ' GP * LAGRANG SOLUTION:',FOPLG C ELSE FOPLG = D0 END IF CALL PRP1AVE(LBESRT(KSYMOP,IOP),AVEVAL,CMO, * WRK(KTUDV),WRK(KWRK1),LWRK1,IPRRSP) HYPFIN = AVEVAL + FOPLG WRITE(LUPRI,'(/2A,3(A,F10.6))') * 'TRIPLET OPERATOR: "',LBESRT(KSYMOP,IOP), * '" LAGRANGIAN:',FOPLG,' AVERAGE:',AVEVAL,' TOTAL:',HYPFIN C C FC and SD values arrays formation C IF (LBESRT(KSYMOP,IOP)(1:3) .EQ. 'FC ') THEN IFCIND=IFCIND+1 HYPVAL(IFCIND,1)=AVEVAL HYPVAL(IFCIND,2)=FOPLG END IF IF (LBESRT(KSYMOP,IOP)(1:3) .EQ. 'SD ') THEN LABEL=LBESRT(KSYMOP,IOP) READ (LABEL,'(3X,I3)') INDSD1 CTMP=LABEL(7:8) indsd2 = -9898 IF (CTMP .EQ. ' x') INDSD2=1 IF (CTMP .EQ. ' y') INDSD2=2 IF (CTMP .EQ. ' z') INDSD2=3 if (indsd2 .eq. -9898) then write(lupri,*) 'hjaaj label=',label, indsd1 write(lupri,*) 'hjaaj no x or y or z in label(7:8) !' call quit('hjaaj: problem in lagrang.F') end if SDVAL(SDIND(2,INDSD1),INDSD2,SDIND(1,INDSD1))=HYPFIN END IF C 100 CONTINUE TRPLET = TRPSAVE END IF DO 300 ISYM = 2,NSYM DO 350 IOP = 1,NESRT(ISYM) WRITE(LUPRI,'(/3A/A,I3)') 'TRIPLET OPERATOR: "', & LBESRT(ISYM,IOP),'" contribution = 0.0 by symmetry.', & '- Operator is of symmetry no.',ISYM 350 CONTINUE 300 CONTINUE C C Singlet operators: C TRPSAVE= TRPLET TRPLET = .FALSE. C ... so PRP1AVE calculates singlet expectation values /hjaaj march 2003 ... C (if TRPLET true, then inactive density matrix is omitted!) DO 400 IOP = 1,NESRS(KSYMOP) CALL PRP1AVE(LBESRS(KSYMOP,IOP),AVEVAL,CMO, * UDV,WRK(KWRK1),LWRK1,IPRRSP) C CALL PRP1AVE(LABELOP,AVEVAL,CMO, UDV,WRK,LWRK,IPRINT) WRITE(LUPRI,'(/3A,F10.6)') * 'SINGLET OPERATOR: "',LBESRS(KSYMOP,IOP),'" AVERAGE:',AVEVAL 400 CONTINUE C DO 500 ISYM = 2,NSYM DO 550 IOP = 1,NESRS(ISYM) WRITE(LUPRI,'(/3A/A,I3)') 'SINGLET OPERATOR: "', & LBESRS(ISYM,IOP),'" contribution = 0.0 by symmetry.', & '- Operator is of symmetry no.',ISYM 550 CONTINUE 500 CONTINUE TRPLET = TRPSAVE C C Isotropic hyperfine coupling printing C IF (FCFLG) THEN CALL AROUND('Isotropic Hyperfine Coupling') WRITE(LUPRI,'(A)') ' Atom Mass G-val ' * //' A, Mhz A, G ' WRITE(LUPRI,'(A)') ' -----------------------------' * //'--------------------------' DO 600 ITMP = 1, ESRNUC DO 610 ISONM = 1, ISODAT(1,NUCINF(ITMP)) IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'GVAL') XISMAS=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'MASS') HYPFIN=((HYPVAL(ITMP,1)+HYPVAL(ITMP,2)) * *XATGV*XTHZ*1.0D-6*ALPHA2) * /(2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) WRITE(LUPRI,'(A,F6.2,A,F8.5,A,F9.4,A,F9.4,A)') * ' * '//NAMN(NUCINF(ITMP))(1:4)//' * ', XISMAS, * ' * ', XATGV,' * ',HYPFIN,' * ',HYPFIN/2.8025D0, * ' *' 610 CONTINUE 600 CONTINUE WRITE(LUPRI,'(A)') ' -----------------------------' * //'--------------------------' WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' * //' symmetry unique atoms!' CALL AROUND('R-U contributions to Isotropic hyperfine' * //' coupling') WRITE(LUPRI,'(A)') ' Atom G-val Average, MHz ' * //' Response, MHz Total, MHz ' WRITE(LUPRI,'(A)') ' ------------------------------' * //'---------------------------------' DO 601 ITMP = 1, ESRNUC DO 611 ISONM = 1, ISODAT(1,NUCINF(ITMP)) IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS,'GVAL') HYPFIN=(XATGV*XTHZ*1.0D-6*ALPHA2) * /(2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) WRITE(LUPRI,'(A,F8.5,A,F11.4,A,F11.4,A,F11.4,A)') * ' * '//NAMN(NUCINF(ITMP))(1:4)//' * ', XATGV, * ' * ',HYPFIN*HYPVAL(ITMP,1),' * ', * HYPFIN*HYPVAL(ITMP,2),' * ', * HYPFIN*(HYPVAL(ITMP,1)+HYPVAL(ITMP,2)),' * ' 611 CONTINUE 601 CONTINUE WRITE(LUPRI,'(A)') ' ------------------------------' * //'---------------------------------' WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' * //' symmetry unique atoms!' ENDIF C C Anisotropic hyperfine coupling printing C IF (SDFLG) THEN CALL AROUND('Anisotropic Hyperfine Coupling') DO 700 ITMP=1,ESRNUC DO 710 ISONM=1,ISODAT(1,NUCINF(ITMP)) IATIS=ISODAT(ISONM+1,NUCINF(ITMP)) XATGV=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS, * 'GVAL') XISMAS=DISOTP(INT(CHARGE(NUCINF(ITMP))),IATIS, * 'MASS') HYPFIN=(XATGV*XTHZ*1.0D-6*ALPHA2)/ * (2.0D0*XPRTMAS*DEGEN(ITMP)*REFSPIN) WRITE(LUPRI,'(/A,F6.2,A,F8.5,A)') * ' A Tensor components for '// * NAMN(NUCINF(ITMP))(1:4)//' ( Mass =', * XISMAS, ' G-val =', XATGV,') in MHz: ' WRITE(LUPRI,'(A)') ' =================================' * //'==================================' WRITE(LUPRI,'(/A)') ' ' WRITE(LUPRI,'(A)') ' SX SY SZ ' * //' ' DO 720 ICRD=1,3 WRITE(LUPRI,'(A,F11.4,A,F11.4,A,F11.4,A)') * ' I'//CHRXYZ(ICRD)//' * ', * HYPFIN*SDVAL(ICRD,1,NUCINF(ITMP)), ' * ', * HYPFIN*SDVAL(ICRD,2,NUCINF(ITMP)), ' * ', * HYPFIN*SDVAL(ICRD,3,NUCINF(ITMP)), ' * ' 720 CONTINUE 710 CONTINUE 700 CONTINUE WRITE(LUPRI,'(/A)') ' *** NOTE: Results printed only for' * //' symmetry unique atoms!' ENDIF C C *** end of RSPESR -- CALL QEXIT('RSPESR') RETURN END C /* Deck solgdt */ SUBROUTINE SOLGDT(CREF,CMO,INDXCI,DV,G,ESOLT,ERLM,WRK, & LFREE,NHCREF,KREFSY) C C Based on SOLGRD: C Copyright 29-Nov-1986 Hans Joergen Aa. Jensen C C Purpose: calculate MCSCF energy and gradient contribution C from a surrounding medium, cavity radius = Rsol C and dielectric constant = EPsol. C C Output: C G MCSCF gradient with solvation contribution added C ESOLT total solvation energy C ERLM(lm,1) contains Esol(l,m) contribution to ESOLT C ERLM(lm,2) contains Tsol(l,m) C #include "implicit.h" #include "dummy.h" C DIMENSION CREF(*), CMO(*), INDXCI(*) DIMENSION DV(*), G(*), ERLM(NLMSOL,2), WRK(*) PARAMETER ( D0 = 0.0D0, D2 = 2.0D0 ) #include "thrzer.h" #include "iratdef.h" #include "priunit.h" #include "infrsp.h" C C Used from common blocks: C INFINP: NLMSOL, LSOLMX, EPSOL, RSOL(3) C INFVAR: NCONF, NWOPT, NVAR, NVARH C INFORB: NNASHX, NNBASX, NNORBX, etc. C INFIND: IROW(*) C INFTAP: LUSOL, LUIT2 C INFPRI: IPRSOL C #include "maxash.h" #include "maxorb.h" #include "infinp.h" #include "infvar.h" #include "inforb.h" #include "infind.h" #include "inftap.h" #include "infpri.h" C LOGICAL FIRST PARAMETER (MXLMAX = 50) DIMENSION ISYRLM(2*MXLMAX+1) CHARACTER*8 STAR8, SOLVDI, EODATA SAVE FIRST DATA FIRST/.TRUE./, STAR8/'********'/ DATA SOLVDI/'SOLVDIAG'/, EODATA/'EODATA '/ C C Statement functions; C define automatic arrays (dynamic core allocation) C FLVEC(LM) = WRK(LM) FLINR(LM) = WRK(KFLINR-1+LM) TLMSI(LM) = WRK(KTLMSI-1+LM) C CALL QENTER('SOLGDT') C IF (LSOLMX .GT. MXLMAX) THEN WRITE (LUERR,*) 'ERROR SOLGDT, increase MXLMAX parameter' WRITE (LUERR,*) ' LSOLMX =',LSOLMX WRITE (LUERR,*) ' MXLMAX =',MXLMAX CALL QUIT('ERROR SOLGDT, increase MXLMAX parameter') END IF C C Core allocation C FLVEC f(l) factors in solvent energy expression C DIASH diagonal of solvent contribution to Hessian C GRDLM TELM gradient for current l,m value in the l,m loop C UCMO CMO unpacked (i.e. no symmetry blocking) C RLMAC active-active subblock of RLM C RLM R(l,m) integrals for current l,m value in l,m loop C C If (INERSF) C (i.e. If (inertial polarization contribution to final state)) C FLINR f(l) factors for inertial pol. contrib. C TLMSI T(lm) values for initial state C end if C KFLVEC = 1 C ... NOTE: KFLVEC = 1 assumed in FLVEC(LM) definition above. IF (INERSF) THEN KFLINR = KFLVEC + NLMSOL KTLMSI = KFLINR + NLMSOL KDIASH = KTLMSI + NLMSOL ELSE KFLINR = 1 KTLMSI = 1 KDIASH = KFLVEC + NLMSOL END IF KGRDLM = KDIASH + NVAR KUCMO = KGRDLM + NVARH KRLMAC = KUCMO + NORBT*NBAST KRLM = KRLMAC + NNASHX KW10 = KRLM + NNORBX C 1.1 read rlmao in ao basis and transform to rlm in mo basis KRLMAO = KW10 KW20 = KRLMAO + NNBASX C 1.2 diagonal contribution for current l,m value in the l,m loop KDIALM = KW10 KW21 = KDIALM + NVAR LW21 = LFREE - KW21 C 1.3 rest of CSF contribution KW22 = KW10 C KTDV = MAX(KW20,KW21,KW22) KWRK1 = KTDV + NASHT * NASHT LWRK1 = LFREE - KWRK1 IF (LWRK1 .LT. 0) CALL ERRWRK('SOLGDT',-KWRK1,LFREE) C IF (IPRSOL .GE. 130) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' SOLGDT - gtot (input) - non-zero elements', * ' NCONF, NWOPT =',NCONF,NWOPT DO 40 I = 1,NCONF IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,G(I) 40 CONTINUE DO 50 I = NCONF+1,NVAR IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,G(I) 50 CONTINUE END IF C C Calculate f(l) factors C If (INERSF) FLVEC factors describe the optical polarization C and FLINR factors describe the inertial polarization C else FLVEC may describe optical or static polarization C CALL SOLFL(WRK(KFLVEC),EPSOL,RSOL,LSOLMX) IF (INERSF) THEN CALL SOLINR(WRK(KFLVEC),WRK(KFLINR),WRK(KTLMSI)) END IF IF ((IPRSOL .GE. 5 .AND. FIRST) .OR. IPRSOL .GE. 15) THEN IF (.NOT.INERSF) THEN WRITE (LUPRI,'(//A/A)') * ' SOLGDT: l f(l) factor', * ' === =================' ELSE WRITE (LUPRI,'(//A/A)') * ' SOLGDT: l optical f(l) factor inertial f(l) factor', * ' === =================== ====================' END IF DO 140 L = 0,LSOLMX LL = (L+1)*(L+1) FL = FLVEC(LL) IF (INERSF) THEN FLI = FLINR(LL) WRITE (LUPRI,'(I15,F17.10,F21.10)') L, FL, FLI ELSE WRITE (LUPRI,'(I15,F16.10)') L, FL END IF 140 CONTINUE END IF C C Read and check dimension information (if first read) and C nuclear contributions to ERLM (always). C CALL GPOPEN(LUSOL,FNSOL,'OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) REWIND LUSOL CALL MOLLAB('SOLVRLM ',LUSOL,LUERR) IF (FIRST) THEN READ (LUSOL) LMAXSS, LMTOT, NNNBAS NERR = 0 IF (LMAXSS .LT. LSOLMX) THEN NERR = NERR + 1 WRITE (LUPRI,'(//2A,2(/A,I5))') ' SOLGDT ERROR,', * ' insufficient number of integrals on LUSOL', * ' l max from SIRIUS input :',LSOLMX, * ' l max from LUSOL file :',LMAXSS END IF IF ((LMAXSS+1)**2 .NE. LMTOT) THEN NERR = NERR + 1 WRITE (LUPRI,'(//2A,3(/A,I5))') ' SOLGDT ERROR,', * ' LUSOL file info inconsistent', * ' l_max :',LMAXSS, * ' (l_max + 1) ** 2 :',(LMAXSS+1)**2, * ' LMTOT :',LMTOT END IF IF (NNNBAS .NE. NBAST) THEN NERR = NERR + 1 WRITE (LUPRI,'(//2A,3(/A,I5))') ' SOLGDT ERROR,', * ' LUSOL file info inconsistent with SIRIUS input', * ' NBAST - LUSOL :',NNNBAS, * ' NBAST - SIRIUS :',NBAST END IF IF (NERR .GT. 0) THEN CALL QUIT('SOLGDT ERROR: LUSOL file not OK for this calc.') END IF ELSE READ (LUSOL) END IF CALL READT(LUSOL,NLMSOL,ERLM(1,2)) C IF (IPRSOL .GE. 20 .AND. NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' SOLGDT - DV matrix :' CALL OUTPAK(DV,NASHT,1,LUPRI) END IF IF (IPRSOL .GE. 7) THEN WRITE (LUPRI,'(/A/)') * ' l, m, Tn(lm) - the nuclear contributions :' LM = 0 DO 220 L = 0,LSOLMX DO 210 M = -L,L LM = LM + 1 WRITE (LUPRI,'(2I5,F15.10)') L,M,ERLM(LM,2) 210 CONTINUE WRITE (LUPRI,'()') 220 CONTINUE END IF C C Unpack symmetry blocked CMO C Loop over l,m expansion C CALL UPKCMO(CMO,WRK(KUCMO)) IF (IPRSOL .GE. 6) * WRITE (LUPRI, '(//A/)') ' SOLGDT: START LOOP OVER LM' CALL DZERO(WRK(KDIASH),NVAR) LM = 0 DO 520 L = 0,LSOLMX READ (LUSOL) L1,(ISYRLM(M),M=1,2*L+1) IF (L1 .NE. L) THEN WRITE (LUERR,*) 'ERROR SOLGDT: L from LUSOL not as expected' WRITE (LUERR,*) 'L from 520 loop:',L WRITE (LUERR,*) 'L from LUSOL :',L1 CALL QUIT('ERROR SOLGDT: L from LUSOL not as expected') END IF DO 500 M = -L,L LM = LM + 1 IF (IPRSOL .GE. 15) THEN WRITE (LUPRI,'(/A,2I5/A)') ' --- l, m :',L,M, * ' ====================' WRITE (LUPRI,'(A,I2)') ' Symmetry :',ABS(ISYRLM(L+M+1)) END IF IF (ISYRLM(L+M+1) .NE. 1) THEN IF (ABS(ERLM(LM,2)) .GT. 1000.D0*THRZER) THEN WRITE (LUPRI,*) 'ERROR SOLGDT for l,m',L,M WRITE (LUPRI,*) 'Symmetry :',ISYRLM(L+M+1) WRITE (LUPRI,*) 'Tn(l,m) .ne. 0, but =',ERLM(LM,2) CALL QUIT('ERROR SOLGDT: Tn(l,m) not 0 as expected') END IF ERLM(LM,2) = D0 C ... to fix round-off errors in Tn(l,m) calculation IF (ISYRLM(L+M+1) .GT. 1) READ (LUSOL) GO TO 500 END IF C C Read R(l,m) in ao basis and transform to mo basis. C Extract active-active block in RLMAC(1) = WRK(KRLMAC). C CALL READT(LUSOL,NNBASX,WRK(KRLMAO)) CALL UTHU(WRK(KRLMAO),WRK(KRLM),WRK(KUCMO),WRK(KWRK1), & NBAST,NORBT) IF (NASHT .GT. 0) THEN CALL GETAC2(WRK(KRLM),WRK(KRLMAC)) END IF IF (IPRSOL .GE. 15) THEN WRITE (LUPRI,'(/A)') ' Rlm_ao matrix:' CALL OUTPAK(WRK(KRLMAO),NBAST,1,LUPRI) WRITE (LUPRI,'(/A)') ' Rlm_mo matrix:' CALL OUTPAK(WRK(KRLM), NORBT,1,LUPRI) IF (NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' Rlm_ac matrix:' CALL OUTPAK(WRK(KRLMAC),NASHT,1,LUPRI) END IF END IF C C Add electronic contribution TE(l,m) to T(l,m) C KFREE=1 ISPIN1=0 ISPIN2=0 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * INDXCI,WRK(KWRK1),KFREE,LWRK1) C C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX C CALL DSITSP(NASHT,WRK(KTDV),DV) C TELM = SOLELM(DV,WRK(KRLMAC),WRK(KRLM),TELMAC) C C construct again triplet density matrix KFREE =1 ISPIN1=1 ISPIN2=0 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * INDXCI,WRK(KWRK1),KFREE,LWRK1) C C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX C CALL DSITSP(NASHT,WRK(KTDV),DV) C IF (IPRSOL .GE. 6) THEN WRITE (LUPRI,'(A,2I5,/A,3F17.8)') * ' --- l, m :',L,M, * ' Te(lm), Tn(lm), T(lm) :', * TELM,ERLM(LM,2),ERLM(LM,2)-TELM IF (IPRSOL .GE. 10) WRITE (LUPRI,'(A,F17.8)') * ' --- active part of Te(lm) :',TELMAC IF (INERSF) WRITE (LUPRI,'(A,F17.8)') * ' --- inertial T(lm) value :',TLMSI(LM) END IF ERLM(LM,2) = ERLM(LM,2) - TELM IF (ABS(ERLM(LM,2)) .LE. THRZER) THEN ERLM(LM,2) = D0 GO TO 500 END IF C ... test introduced 880109 hjaaj C (the only possible problem is the DO 420 loop, C but I think (w.o. having checked) that this C contribution to the Hessian diagonal also will be C zero if ERLM(LM,2) zero). C C Calculate orbital TE(l,m) gradient contribution C and part of csf contribution. C CALL DZERO(WRK(KGRDLM),NVARH) IF (NCONF .GT. 1) THEN CALL SOLGC(CREF,WRK(KRLMAC),TELMAC,WRK(KGRDLM),INDXCI, & WRK(KWRK1),LWRK1) C CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) END IF IF (NWOPT .GT. 0) THEN CALL SOLGO(D0,DV,WRK(KRLM),WRK(KGRDLM+NCONF)) END IF C C C Obtain DIALM = diagonal TE(l,m) Hessian C = 2 ( - TE(l,m) ) C Add the DIALM contribution and the GRDLM contribution C to solvent Hessian diagonal. C Clf the diagonal hessian is not needed for triplet gradients Clf CALL SOLDIA(TELMAC,WRK(KRLMAC),INDXCI, Clf * WRK(KRLM),DV,WRK(KDIALM),WRK(KW21),LW21) C CALL SOLDIA(TELM,RLMAC,INDXCI,RLM,DV,DIAG,WRK,LFREE) C FAC1 = - D2 * FLVEC(LM) * ERLM(LM,2) IF (INERSF) THEN FAC1 = FAC1 - FLINR(LM) * D2 * TLMSI(LM) END IF FAC2 = D2 * FLVEC(LM) DO 420 I = 0,(NVAR-1) WRK(KDIASH+I) = WRK(KDIASH+I) * + FAC1 * WRK(KDIALM+I) * + FAC2 * WRK(KGRDLM+I) * WRK(KGRDLM+I) 420 CONTINUE C C test orthogonality C IF (IPRSOL .GE. 120) THEN WRITE (LUPRI,'(/A)')' SOLGDT - grdlm, dialm, diash, cref' DO 430 I = 1,NCONF WRITE (LUPRI,'(A,I10,4F12.6)') ' conf #',I, * WRK(KGRDLM-1+I),WRK(KDIALM-1+I),WRK(KDIASH-1+I),CREF(I) 430 CONTINUE END IF TEST = DDOT(NCONF,CREF,1,WRK(KGRDLM),1) IF (ABS(TEST) .GT. 1.D-8) THEN NWARN = NWARN + 1 WRITE (LUPRI,'(/A,I5,/A,1P,D12.4)') * ' SOLGDT WARNING, for LM =',LM, * ' =',TEST END IF C C Add TE(l,m) gradient contribution to MCSCF gradient C g = g - 2 f(l) * T(l,m) * (dTE(l,m)/d(lambda)) C FAC = - D2 * FLVEC(LM) * ERLM(LM,2) IF (INERSF) THEN FAC = FAC - FLINR(LM) * D2 * TLMSI(LM) END IF CALL DAXPY(NVARH,FAC,WRK(KGRDLM),1,G,1) IF (IPRSOL .GE. 140) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' SOLGDT - grdlm, gtot (accum) - non-zero grdlm', * ' NCONF, NWOPT =',NCONF,NWOPT DO 440 I = 1,NCONF IF (WRK(KGRDLM-1+I) .NE. D0) * WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,FAC*WRK(KGRDLM-1+I),G(I) 440 CONTINUE DO 450 I = NCONF+1,NVAR IF (WRK(KGRDLM-1+I) .NE. D0) * WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,FAC*WRK(KGRDLM-1+I),G(I) 450 CONTINUE END IF C 500 CONTINUE 520 CONTINUE C CALL GPCLOSE(LUSOL,'KEEP') C C 500 is end of (l,m) loop. C C IF (IPRSOL .GE. 130) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' SOLGDT - gtot (output) - non-zero elements', * ' NCONF, NWOPT =',NCONF,NWOPT DO 840 I = 1,NCONF IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,G(I) 840 CONTINUE DO 850 I = NCONF+1,NVAR IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,G(I) 850 CONTINUE END IF C C C Calculate ER(l,m) energy contributions and add them up C ESOLT = D0 DO 900 LM = 1,NLMSOL ERLM(LM,1) = FLVEC(LM) * ERLM(LM,2) * ERLM(LM,2) IF (INERSF) THEN ERLM(LM,1) = ERLM(LM,1) * + FLINR(LM) * ERLM(LM,2) * D2 * TLMSI(LM) * - FLINR(LM) * TLMSI(LM) * TLMSI(LM) END IF ESOLT = ESOLT + ERLM(LM,1) 900 CONTINUE C FIRST = .FALSE. CALL QEXIT('SOLGDT') RETURN C end of solgdt. END C /* Deck FCOPER */ SUBROUTINE FCOPER(ATMIND,LABINT) #include "implicit.h" #include "priunit.h" #include "mxcent.h" CHARACTER*8 LABINT INTEGER ATMIND #include "nuclei.h" #include "chrnos.h" C LABINT = 'FC '//NAMN(ATMIND)(1:2) & //CHRNOS(ATMIND/100) & //CHRNOS(MOD(ATMIND,100)/10) & //CHRNOS(MOD(ATMIND,10)) C RETURN END C /* Deck SETDEG */ SUBROUTINE SETDEG #include "implicit.h" #include "priunit.h" #include "maxaqn.h" #include "maxmom.h" #include "mxcent.h" #include "maxorb.h" #include "parnmr.h" #include "nuclei.h" #include "symmet.h" #include "chrnos.h" #include "chrxyz.h" C DO 50 I = 1, ESRNUC DEGEN(I) = 0.0D0 DO 100 IREP=0,MAXREP ISCOR1 = IPTCNT(3*(NUCINF(I) - 1) + 1,IREP,2) IF (ISCOR1 .GT. 0) DEGEN(I)=DEGEN(I)+1.0D0 100 CONTINUE 50 CONTINUE C RETURN END C /* Deck SETISO */ SUBROUTINE SETISO #include "implicit.h" #include "priunit.h" #include "parnmr.h" #include "mxcent.h" #include "nuclei.h" C DO 100 I=1,ESRNUC IF (ISODAT(1,NUCINF(I)) .EQ. 0) THEN ISODAT(1,NUCINF(I))=1 IATIS=1 200 XATGV=DISOTP(INT(CHARGE(NUCINF(I))),IATIS,'GVAL') IF (DABS(XATGV) .LT. 1.0D-5) THEN IATIS=IATIS+1 GO TO 200 ENDIF ISODAT(2,NUCINF(I))=IATIS ENDIF 100 CONTINUE C RETURN END C /* Deck SDOPER */ SUBROUTINE SDOPER(ATMIND,LABINT,NCOOR) #include "implicit.h" #include "priunit.h" #include "maxaqn.h" #include "maxmom.h" #include "mxcent.h" #include "maxorb.h" #include "parnmr.h" INTEGER ATMIND, NCOOR CHARACTER LABINT(9*ATMNUM)*8 #include "nuclei.h" #include "symmet.h" #include "chrnos.h" #include "chrxyz.h" C NCOOR=0 DEGEN(ATMIND) = 0.0D0 DO 100 IREP=0,MAXREP DO 200 ICOOR1=1,3 ISCOR1 = IPTCNT(3*(NUCINF(ATMIND ) - 1) + ICOOR1,IREP,2) IF (ISCOR1 .GT. 0) THEN IF (ICOOR1 .EQ. 1) DEGEN(ATMIND)=DEGEN(ATMIND)+1.0D0 DO 300 ICOOR2 = 1, 3 NCOOR=NCOOR+1 LABINT(NCOOR) = 'SD '//CHRNOS(ISCOR1/100) & //CHRNOS(MOD(ISCOR1,100)/10) & //CHRNOS(MOD(ISCOR1,10)) & //' '//CHRXYZ(-ICOOR2) 300 CONTINUE SDIND(1,ISCOR1)=ATMIND SDIND(2,ISCOR1)=ICOOR1 ENDIF 200 CONTINUE 100 CONTINUE C RETURN END C /* Deck pcmgdt */ SUBROUTINE PCMGDT(CREF,CMO,INDXCI,DV,G,ESOLT,WRK,LFREE, $ NHCREF,KREFSY) C C Based on SOLGDT: C 13-02-2006 Luca Frediani C C Purpose: calculate MCSCF energy and gradient contribution C from a surrounding medium with PCM C C Output: C G MCSCF gradient with solvation contribution added C #include "implicit.h" #include "dummy.h" C DIMENSION CREF(*), CMO(*), INDXCI(*) DIMENSION DV(*), G(*), WRK(*) PARAMETER ( D1=1.0d0, DM1=-1.0d0, D0 = 0.0D0, D2 = 2.0D0 ) #include "thrzer.h" #include "iratdef.h" #include "priunit.h" #include "infrsp.h" #include "mxcent.h" #include "orgcom.h" #include "pcmdef.h" #include "pcmlog.h" #include "pcm.h" C C Used from common blocks: C INFINP: NLMSOL, LSOLMX, EPSOL, RSOL(3) C INFVAR: NCONF, NWOPT, NVAR, NVARH C INFORB: NNASHX, NNBASX, NNORBX, etc. C INFIND: IROW(*) C INFTAP: LUSOL, LUIT2 C INFPRI: IPRSOL C #include "maxash.h" #include "maxorb.h" #include "infinp.h" #include "infvar.h" #include "inforb.h" #include "infind.h" #include "inftap.h" #include "infpri.h" C LOGICAL FNDLAB, EXP1VL, TOFILE, TRIMAT PARAMETER (MXLMAX = 50) CHARACTER*8 STAR8, SOLVDI, EODATA, LABINT(9*MXCENT) DATA STAR8/'********'/ DATA SOLVDI/'SOLVDIAG'/, EODATA/'EODATA '/ C C Statement functions; C define automatic arrays (dynamic core allocation) C CALL QENTER('PCMGDT') C C C Core allocation C DIASH diagonal of solvent contribution to Hessian C GRDLM TELM gradient for current l,m value in the l,m loop C UCMO CMO unpacked (i.e. no symmetry blocking) C KJENAO = 1 KJEN = KJENAO + NNBASX KJ1AO = KJEN + NNORBX KJ1 = KJ1AO + NNBASX*NSYM KJENAC = KJ1 + NNORBX KJ1AC = KJENAC + NNASHX KDIASH = KJ1AC + NNASHX KGRDLM = KDIASH + NVAR KUCMO = KGRDLM + NVARH KJ2GRD = KUCMO + NORBT*NBAST KPOT = KJ2GRD + NVAR KDENC = KPOT + NTS KDENV = KDENC + N2BASX KW10 = KDENV + N2BASX C 1.3 rest of CSF contribution KDIALM = KW10 KTDV = KDIALM + NVAR KW20 = KTDV + NASHT * NASHT C C Allocations for non-equlibrium energy solvation IF (NONEQ) THEN KMPOT = KW20 KQSEGR = KMPOT + NTS * NTS KPOTGR = KQSEGR + NTS KW21 = KPOTGR + NTS ELSE KW21 = KW20 END IF C LW21 = LFREE - KW21 C KWRK1 = KW21 LWRK1 = LFREE - KWRK1 IF (LWRK1 .LT. 0) CALL ERRWRK('PCMGDT',-KWRK1,LFREE) C KFREE=1 ISPIN1=0 ISPIN2=0 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * INDXCI,WRK(KWRK1),KFREE,LWRK1) C C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX C CALL DSITSP(NASHT,WRK(KTDV),DV) C clf IF (IPRPCM .GE. 130) THEN IF (.true.) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' --- PCMGDT - gtot (input) - non-zero elements', * ' NCONF, NWOPT =',NCONF,NWOPT DO 40 I = 1,NCONF IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,G(I) 40 CONTINUE DO 50 I = NCONF+1,NVAR IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,G(I) 50 CONTINUE END IF IF (IPRSOL .GE. 20 .AND. NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' SOLGDT - DV matrix :' CALL OUTPAK(DV,NASHT,1,LUPRI) END IF C C Unpack symmetry blocked CMO C Loop over l,m expansion C CALL UPKCMO(CMO,WRK(KUCMO)) CALL DZERO(WRK(KDIASH),NVAR) CALL DZERO(WRK(KDIALM),NVAR) CALL DZERO(WRK(KJEN),NNORBX) C C Read JEN = J(en) + J(ne) in ao basis and transform to mo basis. C Extract active-active block in WRK(KJENAC). C LUPBKP = LUPROP IF (LUPROP .LT. 0) CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN', & 'SEQUENTIAL','UNFORMATTED',IDUMMY,.FALSE.) CLF IF (.NOT. (FNDLAB('NE-PCMIN',LUPROP))) THEN IF (.TRUE.) THEN CALL PCMJMAT(WRK(KJENAO),NNBASX,WRK(KWRK1),LWRK1) END IF REWIND (LUPROP) CALL REAPCM('NE-PCMIN','PCMGRD ',LUPROP,WRK(KJENAO),NNBASX) CALL UTHU(WRK(KJENAO),WRK(KJEN),WRK(KUCMO),WRK(KWRK1), & NBAST,NORBT) IF (NASHT .GT. 0) THEN CALL GETAC2(WRK(KJEN),WRK(KJENAC)) END IF IF (IPRPCM .GE. 15) THEN WRITE (LUPRI,'(/A)') ' JEN_ao matrix:' CALL OUTPAK(WRK(KJENAO),NBAST,1,LUPRI) WRITE (LUPRI,'(/A)') ' JEN_mo matrix:' CALL OUTPAK(WRK(KJEN), NORBT,1,LUPRI) IF (NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' JEN_ac matrix:' CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) END IF END IF C C Expextation value of JEN (=PB) C TJEN = SOLELM(DV,WRK(KJENAC),WRK(KJEN),TJENAC) IF (IPRPCM .GE. 6) THEN WRITE (LUPRI,'(A,F17.8)') * ' --- JEN expectation value(=PB) :',TJEN WRITE (LUPRI,'(A,F17.8)') * ' --- active part of JEN(=PB) :',TJENAC END IF Cbm PB=-TJEN Cbm WRITE(LUPRI,*)'PB =',PB C C Read J2 in ao basis and transform to mo basis. C Extract active-active block in WRK(KJ2AC). C CALL DZERO(WRK(KDENC),N2BASX) CALL DZERO(WRK(KDENV),N2BASX) CALL FCKDEN((NISHT.GT.0),(NASHT.GT.0),WRK(KDENC),WRK(KDENV), & CMO,DV,WRK(KWRK1),LW21) CALL DAXPY(N2BASX,1.0D0,WRK(KDENV),1,WRK(KDENC),1) CALL DZERO(WRK(KDENV),N2BASX) CALL DGEFSP(NBAST,WRK(KDENC),WRK(KDENV)) ckr CALL DZERO(WRK(KDENC),N2BASX) ckr CALL PKSYM1(WRK(KDENV),WRK(KDENC),NBAS,NSYM,1) EXP1VL = .TRUE. TOFILE = .FALSE. CALL J1INT(WRK(KPOT),EXP1VL,WRK(KDENV),1,TOFILE,'NPETES ', & 1,WRK(KWRK1),LW21) CALL DZERO(QSE,MXTS) CALL V2Q(WRK(KWRK1),WRK(KPOT),QSE,QET,.FALSE.) CALL GPCLOSE(LUPCMD,'KEEP') C CALL DSCAL(NTS,-1.0D0,QSE,1) C C Read J1 (=ELECTROSTATIC POTENTIAL) in ao basis and transform to mo C basis. Extract active-active block in WRK(KJ1AC). C XI = DIPORG(1) YI = DIPORG(2) ZI = DIPORG(3) Ckr Ckr should be parallelized, but a lot of information to send Ckr However, in general not used for SCF and DFT optimizations Ckr DO I = 1 , NTSIRR L = 1 NCOMP = NSYM DIPORG(1) = XTSCOR(I) DIPORG(2) = YTSCOR(I) DIPORG(3) = ZTSCOR(I) EXP1VL = .FALSE. TOFILE = .FALSE. KPATOM = 0 TRIMAT = .TRUE. CALL GET1IN(WRK(KJ1AO),'NPETES ',NCOMP,WRK(KWRK1),LW21,LABINT, & INTREP,INTADR,L,TOFILE,KPATOM,TRIMAT,DUMMY,EXP1VL, & DUMMY,IPRPCM) CALL UTHU(WRK(KJ1AO),WRK(KJ1),WRK(KUCMO),WRK(KWRK1), & NBAST,NORBT) IF (NASHT .GT. 0) THEN CALL GETAC2(WRK(KJ1),WRK(KJ1AC)) END IF IF (NONEQ) THEN Clf here we will need the singlet density, I beleive..... WRK(KPOT + I - 1) = SOLELM(DV,WRK(KJ1AC),WRK(KJ1),TJ1AC) END IF IF (IPRPCM .GE. 15) THEN WRITE (LUPRI,'(/A)') ' J1_ao matrix:' CALL OUTPAK(WRK(KJ1AO),NBAST,1,LUPRI) WRITE (LUPRI,'(/A)') ' J1_mo matrix:' CALL OUTPAK(WRK(KJ1), NORBT,1,LUPRI) IF (NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' J1_ac matrix:' CALL OUTPAK(WRK(KJ1AC),NASHT,1,LUPRI) END IF END IF CALL DAXPY(NNORBX,-QSE(I),WRK(KJ1),1,WRK(KJEN),1) C C Due to the computational cost of building the CI gradient, and the lack C of importance on convergence rates, we skip the construction of the C solvent contribution to the diagonal Hessian. K.Ruud, Oct.-01 C C CALL DZERO(WRK(KGRDLM),NVARH) C IF (NCONF .GT. 1) THEN C CALL SOLGC(CREF,WRK(KJ1AC),TJ1AC,WRK(KGRDLM),INDXCI, C & WRK(KWRK1),LWRK1) CC CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) CC END IF C IF (NWOPT .GT. 0) THEN C CALL SOLGO(DCVAL,DV,WRK(KJ1),WRK(KGRDLM+NCONF)) C END IF C READ(LUGRDQ)WRK(KJ2GRD) C DO J = NCONF, NVAR - 1 C WRK(KDIASH+J) = WRK(KDIASH+J) - WRK(KGRDLM+J)*WRK(KJ2GRD+J) C ENDDO ENDDO Ckr Ckr End of parallelization loop Ckr C C Expextation value of J + X(0) = PB + PX C IF (NASHT .GT. 0) THEN CALL GETAC2(WRK(KJEN),WRK(KJENAC)) END IF TJEN = SOLELM(DV,WRK(KJENAC),WRK(KJEN),TJENAC) IF (IPRPCM .GE. 6) THEN CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) WRITE (LUPRI,'(A,F17.8)') * ' --- TJEN expectation value :',TJEN WRITE (LUPRI,'(A,F17.8)') * ' --- active part of TJEN :',TJENAC END IF C KFREE=1 ISPIN1=1 ISPIN2=0 CALL RSPDM(KREFSY,KREFSY,NHCREF,NHCREF,CREF,CREF, * WRK(KTDV),DUMMY, ISPIN1,ISPIN2,.FALSE.,.TRUE., * INDXCI,WRK(KWRK1),KFREE,LWRK1) Clf IF (.true.) THEN WRITE (LUPRI,'(/A)') ' JEN_ao matrix:' CALL OUTPAK(WRK(KJENAO),NBAST,1,LUPRI) WRITE (LUPRI,'(/A)') ' JEN_mo matrix:' CALL OUTPAK(WRK(KJEN), NORBT,1,LUPRI) IF (NASHT .GT. 0) THEN WRITE (LUPRI,'(/A)') ' JEN_ac matrix:' CALL OUTPAK(WRK(KJENAC),NASHT,1,LUPRI) END IF END IF C C TRIANGULAR PACKING OF ONE ELECTRON DENSITY MATRIX C CALL DSITSP(NASHT,WRK(KTDV),DV) CALL DZERO(WRK(KGRDLM),NVARH) IF (NCONF .GT. 1) THEN CALL SOLGC(CREF,WRK(KJENAC),TJENAC,WRK(KGRDLM),INDXCI, & WRK(KWRK1),LWRK1) C CALL SOLGC(CREF,RLMAC,TELMAC,GLMCI,INDXCI,WRK,LWRK) END IF IF (NWOPT .GT. 0) THEN CALL SOLGO(D0,DV,WRK(KJEN),WRK(KGRDLM+NCONF)) END IF C C C Obtain DIALM = diagonal TE(l,m) Hessian C = 2 ( - TE(l,m) ) C Add the DIALM contribution and the GRDLM contribution C to solvent Hessian diagonal. C C CALL SOLDIA(TJENAC,WRK(KJENAC),INDXCI, C * WRK(KJEN),DV,WRK(KDIALM),WRK(KW21),LW21) C CALL SOLDIA(TELM,RLMAC,INDXCI,RLM,DV,DIAG,WRK,LFREE) C DO 420 I = 0,(NVAR-1) WRK(KDIASH+I) = WRK(KDIASH+I) * - WRK(KDIALM+I) 420 CONTINUE C C test orthogonality C IF (IPRPCM .GE. 120) THEN WRITE (LUPRI,'(/A)')' --- PCMGRD - grdj1, grdj2, dialm, '// & 'diash, cref' DO 430 I = 1,NCONF WRITE (LUPRI,'(A,I10,3F10.6)') ' conf #',I, * WRK(KDIALM-1+I), * WRK(KDIASH-1+I),CREF(I) 430 CONTINUE END IF C TEST = DDOT(NCONF,CREF,1,WRK(KGRDLM),1) IF (ABS(TEST) .GT. 1.D-8) THEN NWARN = NWARN + 1 WRITE (LUPRI,'(/A,/A,1P,D12.4)') * ' --- PCMGRD WARNING, for B', * ' =',TEST END IF C C TEST = DDOT(NCONF,CREF,1,WRK(KGRDJ1),1) IF (ABS(TEST) .GT. 1.D-8) THEN NWARN = NWARN + 1 WRITE (LUPRI,'(/A,/A,1P,D12.4)') * ' --- PCMGRD WARNING, for J1 ', * ' =',TEST END IF C TEST = DDOT(NCONF,CREF,1,WRK(KGRDJ2),1) IF (ABS(TEST) .GT. 1.D-8) THEN NWARN = NWARN + 1 WRITE (LUPRI,'(/A,/A,1P,D12.4)') * ' --- PCMGRD WARNING, for J2', * ' =',TEST END IF C C Add PCM gradient contribution to MCSCF gradient C CALL DAXPY(NVARH,DM1,WRK(KGRDLM),1,G,1) clf IF (IPRPCM .GE. 140) THEN IF (.true.) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' --- PCMGRD - grdB, gtot (accum) - non-zero grdlm', * ' NCONF, NWOPT =',NCONF,NWOPT DO 440 I = 1,NCONF IF (WRK(KGRDLM-1+I) .NE. D0) * WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,WRK(KGRDLM-1+I),G(I) 440 CONTINUE DO 450 I = NCONF+1,NVAR IF (WRK(KGRDLM-1+I) .NE. D0) * WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,WRK(KGRDLM-1+I),G(I) 450 CONTINUE END IF C C C IF (IPRPCM .GE. 130) THEN WRITE (LUPRI,'(/A/A,2I10)') * ' --- PCMGRD - gtot (output) - non-zero elements', * ' NCONF, NWOPT =',NCONF,NWOPT DO 840 I = 1,NCONF IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' conf #',I,G(I) 840 CONTINUE DO 850 I = NCONF+1,NVAR IF (G(I) .NE. D0) WRITE (LUPRI,'(A,I10,3F15.10)') * ' orb #',I,G(I) 850 CONTINUE END IF IF (LUPBKP .LT. 0) CALL GPCLOSE(LUPROP,'KEEP') CALL QEXIT('PCMGDT') RETURN C end of pcmgdt. END ! --- end of DALTON/rsp/lagrang.F ---