C -*- Fortran -*-
C DELCD DEBUG
C
C..DELCN FORMAT MICHEL 12/2003
C..DELCI CONSTRUCTION OF BIELE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ORTHO - programs for ab-initio calculations in localized orbitals
C Copyright (C) 2008 Peter Reinhardt (Univ. Paris VI, France)
C
C This program is free software: you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation, either version 3 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program. If not, see .
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C..FILE 'flow.h'
C..FILE 'common_syst.h'
C..FILE 'common_bas.h'
C..FILE 'common_intu.h'
C
C..FILE 'common_boys.h'
C
C..FILE 'common_energy.h'
C..FILE 'common_units.h'
C..FILE 'common_vect.h'
C
C..FILE 'common_mezey.h'
C..FILE 'common_two.h'
C..FILE 'common_oda.h'
C
PROGRAM ORS
C
C one single molecule
C
C integral files are KINETIC, OVERLAP, HAMILTO, AOTWO
C
C we read the starting vector
C
C Diagonalization of the CI matrix is done with the Davidson
C algorithm
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C P.Reinhardt, 3/2001
C 3/2003 Boys localization
C 2/2005 Optimal damping algorithm
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCLUDE 'param.h'
PARAMETER (NNLMX=NLMAX+NLMAX-1)
C
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
COMMON /UNITS/ IUNITR,IUNITX,IUNITO
C.. INCLUDE 'common_units.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
C we have the "old" best density and the new density
C the density is PODA
COMMON /ODACOM/ PODA(NBASM,NBASM),XIODA(NBASM,NBASM),EMODA,EKODA
$ ,EHODA,EMPODA,LODAU
LOGICAL LODAU
C.. INCLUDE 'common_oda.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
C
CHARACTER*3 CLSTRG(0:3,7),CBSTR(NBASM)
CHARACTER*8 PNAME
CHARACTER*1 BBBB
PNAME='O R S '
INCLUDE 'compiler_stamp'
X=CPTIME(3)
CALL DATING(PNAME,1)
C
CLSTRG(0,1)='S '
CLSTRG(1,1)='PX '
CLSTRG(1,2)='PY '
CLSTRG(1,3)='PZ '
CLSTRG(2,1)='D-2'
CLSTRG(2,2)='D-1'
CLSTRG(2,3)='D 0'
CLSTRG(2,4)='D 1'
CLSTRG(2,5)='D 2'
CLSTRG(3,1)='F-3'
CLSTRG(3,2)='F-2'
CLSTRG(3,3)='F-1'
CLSTRG(3,4)='F 0'
CLSTRG(3,5)='F 1'
CLSTRG(3,6)='F 2'
CLSTRG(3,7)='F 3'
C
IUNITR=26
IUNITO=15
IUNITX=31
C
ANTOAU=0.529177247D0
EPS=1.D-14
EPSE=1.D-11
NITMAX=50
C
WRITE(6,*)
WRITE(6,*) ' O R S - Version for molecules'
WRITE(6,*)
WRITE(6,*) ' P.Reinhardt, Paris 3/2001 '
WRITE(6,*)
WRITE(6,*)
WRITE(6,*)
C
C read flow options from File INPUT.ORS
C
CALL RDORS
C
C print result of options scan
C
IF (LCAN) THEN
WRITE(6,*) ' CANONICAL SCF PROCEDURE '
IF (LCOREH) THEN
WRITE(6,*) ' Starting iterations with the core Hamiltonian'
ELSE
WRITE(6,*)
$ ' Using the complete starting vector for the iterations'
END IF
ELSE
LCOREH=.FALSE.
WRITE(6,*) ' PRODUCING LOCALIZED ORBITALS '
IF (LCICOM) THEN
WRITE(6,*)
WRITE(6,*) ' constructing the complete SCI matrix'
$ ,' with all bielectronic integrals'
WRITE(6,*)
END IF
END IF
IF (LECONV) THEN
WRITE(6,*) ' WE CONVERGE ON ENERGY'
ELSE
WRITE(6,*) ' WE CONVERGE ON SFN'
END IF
WRITE(6,*) ' WE TRY AT MOST ',NITMAX
$ ,' ITERATIONS IN THE SCF LOOP '
C
C the a posteriori localization options
C
IF (LPIPM.OR.LBOYS) THEN
IF (LPIPM) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE PIPEK-MEZEY LOCALIZED ORBITALS '
WRITE(6,*)
END IF
IF (LBOYS) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE BOYS LOCALIZED ORBITALS '
WRITE(6,*)
END IF
IF (NFRZ.EQ.0) THEN
WRITE(6,*) ' NO ORBITAL FREEZING IN LOCALIZATION PROCEDURE'
ELSE
IF (NFRZ.EQ.1) THEN
WRITE(6,*) NFRZ,' ORBITAL FROZEN IN LOCALIZATION PROCEDURE'
ELSE
WRITE(6,*) ' ',NFRZ
$ ,' ORBITALS FROZEN IN LOCALIZATION PROCEDURE'
WRITE(6,*) ' ORBITALS No ',(IFRZ(J),J=1,NFRZ)
END IF
END IF
WRITE(6,*)
WRITE(6,*) ' for the virtual space '
IF (LREMO) THEN
WRITE(6,*) ' we remove atomic orbitals via a predefined list '
WRITE(6,*)
ELSE
WRITE(6,*)
$ ' we remove atomic orbitals via the overlap criterion'
WRITE(6,*)
END IF
ELSE
NFRZ=0
END IF
IF (LAMBDA) THEN
WRITE(6,*)
WRITE(6,*) ' SCALED BIELECTRONIC INTEGRALS: LAMBDA = ',XLAMB
WRITE(6,*)
END IF
IF (LEXTER) THEN
WRITE(6,*)
WRITE(6,*) ' External construction of the Fock matrix '
WRITE(6,*)
OPEN(UNIT=45,FILE='molpro.out',STATUS='UNKNOWN')
CLOSE(45,STATUS='DELETE')
WRITE(6,*)
WRITE(6,*)
$ ' weight of exact exchange in the density functional : '
$ ,XFAC
WRITE(6,*)
END IF
IF (LODA) THEN
WRITE(6,*) ' The optimal damping algorithm of E.Cances '
END IF
IF (LWFN) THEN
WRITE(6,*)
write(6,*) ' NOT YET IMPLEMENTED, WORK HARDER '
WRITE(6,*) ' Writing the WFN file in GAUSSIAN style to '
WRITE(6,*)
STOP
END IF
IF (LMULTP) THEN
WRITE(6,*)
WRITE(6,*) ' we use additional one-electron integrals from the '
WRITE(6,*) ' file '
WRITE(6,*)
END IF
C
IF (L6D) THEN
WRITE(6,*)
WRITE(6,*) ' we calculate in cartesina Gaussians '
WRITE(6,*) ' ATTENTION : reduced functionality '
WRITE(6,*) ' no localization '
WRITE(6,*)
END IF
C
IF (LICSCF) THEN
WRITE(6,*)
WRITE(6,*) ' including Davidson'' ICSCF (JCP 57 (1972) 1999) '
WRITE(6,*) ' we need canonical orbitals for that '
WRITE(6,*)
END IF
WRITE(6,*)
WRITE(6,*) ' SHIFT OF THE ATOMIC CO-ORDINATES (in a.u.)'
WRITE(6,'(3F16.8)') (ORIGIN(I),I=1,3)
WRITE(6,*)
C read information on system
C
C we read the basis correctly and in all detail, once and for all
C
OPEN(UNIT=IUNITR,FILE='SYSTEM.ORTHO',STATUS='OLD',
- FORM='FORMATTED',ERR=901)
READ(IUNITR,*) NATOM
WRITE(6,*)
C
C MOLCAS, DALTON: first all s, then all p, then all d ...
C so read (sort has been done by GENINPUT) given basis set
C
C well, we have to recalculate the nuclear repulsion of the molecule
C
NSHL=0
NBAS=1
LMAX=0
IPRIM=0
DO IAT=1,NATOM
READ(IUNITR,*) NZ(IAT),NSH(IAT),(POS(J,IAT),J=1,3)
ISHST(IAT)=NSHL+1
IBFST(IAT)=NBAS
DO ICOORD=1,3
POS(ICOORD,IAT)=POS(ICOORD,IAT)-ORIGIN(ICOORD)
END DO
DO ISH=1,NSH(IAT)
NSHL=NSHL+1
READ(IUNITR,*) ITYPE,NPRIM(NSHL)
IL(NSHL)=ITYPE
LMAX=MAX(LMAX,ITYPE+1)
IF (L6D) THEN
IF (ITYPE.EQ.0) THEN
NBAS=NBAS+1
ELSE IF (ITYPE.EQ.1) THEN
NBAS=NBAS+3
ELSE IF (ITYPE.EQ.2) THEN
NBAS=NBAS+6
ELSE IF (IYTPE.EQ.3) THEN
NBAS=NBAS+10
ELSE IF (ITYPE.EQ.4) THEN
NBAS=NBAS+15
ELSE
WRITE(6,*) ' no functions beyond g available, sorry '
STOP
END IF
ELSE
DO III=1,ITYPE+ITYPE+1
C store the atomic centers of the basis functions (needed for Pipek
c -Mezey)
IATZ(NBAS)=IAT
NBAS=NBAS+1
END DO
END IF
DO III=1,NPRIM(NSHL)
IPRIM=IPRIM+1
READ(IUNITR,*) EXX(IPRIM),COEFF(IPRIM)
END DO
END DO
NFIN(IAT)=NBAS-1
END DO
NBAS=NBAS-1
CLOSE(IUNITR)
C
C fill IPERM for the external construction
C MOLPRO orders the d orbitals as 0 -2 1 2 -1
C we have them as -2 -1 0 1 2
C we may have to reorder the atoms
IF (.NOT.LREORD) THEN
DO IAT=1,NATOM
IATSRT(IAT)=IAT
END DO
END IF
DO IALPH=1,NBAS
ISIG(IALPH)=1
END DO
C
C we should be able to count the basis
C -each atom has NSH(IAT) Shells, starting at ISHST(IAT), of type
c IL(ISHL)
C
C we have to assemble an array indicating where the MOLPRO basis
C functions go to
C
C Atoms:
C MOLPRO ORTHO
C 1 -> IATSRT(1)
C ... ...
C NATOM -> IATSRT(NATOM)
C
IF (LEXTER) THEN
IPOS=0
DO IIAT=1,NATOM
IAT=IATSRT(IIAT)
WRITE(6,*) 'MOLPRO Atom No',IIAT,' goes to ORTHO Atom No'
$ ,IAT
IBF=IBFST(IAT)
WRITE(6,*) ' we start at AO No ',IBF
DO ISHL=ISHST(IAT),ISHST(IAT)+NSH(IAT)-1
ITYPE=IL(ISHL)
IF (ITYPE.EQ.0) THEN
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(0,1)
IPERM(IPOS)=IBF
IBF=IBF+1
ELSE IF (ITYPE.EQ.1) THEN
DO IM=1,3
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(1,IM)
IPERM(IPOS)=IBF
IBF=IBF+1
END DO
ELSE IF (ITYPE.EQ.2) THEN
C MOLPRO orders the d oritals as 0 -2 1 2 -1
C 3 1 4 5 2
C we have them as -2 -1 0 1 2
C
IPOS=IPOS+1
IBF=IBF-1
CBSTR(IPOS)=CLSTRG(2,3)
IPERM(IPOS)=IBF+3
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(2,1)
IPERM(IPOS)=IBF+1
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(2,4)
IPERM(IPOS)=IBF+4
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(2,5)
IPERM(IPOS)=IBF+5
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(2,2)
IPERM(IPOS)=IBF+2
IBF=IBF+6
ELSE IF (ITYPE.EQ.3) THEN
C MOLPRO orders the f oritals as 1 -1 0 3 -2 -3 2
C 5 3 4 7 2 1 6
C we have them as -3 -2 -1 0 1 2 3
IBF=IBF-1
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(3,5)
IPERM(IPOS)=IBF+5
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(3,3)
IPERM(IPOS)=IBF+3
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(3,4)
IPERM(IPOS)=IBF+4
IPOS=IPOS+1
ISIG(IPOS)=-1
CBSTR(IPOS)=CLSTRG(3,7)
IPERM(IPOS)=IBF+7
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(3,2)
IPERM(IPOS)=IBF+2
IPOS=IPOS+1
CBSTR(IPOS)=CLSTRG(3,1)
IPERM(IPOS)=IBF+1
IPOS=IPOS+1
ISIG(IPOS)=-1
CBSTR(IPOS)=CLSTRG(3,6)
IPERM(IPOS)=IBF+6
IBF=IBF+8
ELSE
WRITE(6,*) ' g functions and beyond not yet available '
STOP 'work harder '
END IF
END DO
END DO
WRITE(6,*)
DO I=1,NBAS
WRITE(6,9226) I,CBSTR(I),ISIG(IPERM(I))*IPERM(I)
9226 FORMAT(' MOLPRO AO No ',I5,' (',A3,') goes to ORTHO AO No ',I4)
END DO
END IF
C
C again for Pipek-Mezey two auxiliary arrays
C
IAT=1
NABAS(IAT)=0
ISTRT(IAT)=1
DO I=1,NBAS
IF (IATZ(I).EQ.IAT) THEN
NABAS(IAT)=NABAS(IAT)+1
ELSE
IAT=IAT+1
ISTRT(IAT)=I
NABAS(IAT)=1
END IF
END DO
DO IAT=1,NATOM
NABAS(IAT)=NABAS(IAT)-1
IF (LPIPM.OR.LBOYS) WRITE(6,*)
$ ' for Pipek-Mezey: ATOM, ISTRT, NABAS ',IAT, ISTRT(IAT),
$ NABAS(IAT)
END DO
C
IF (NFRZ.NE.0) THEN
DO I=1,NBAS
ITMP(I)=0
END DO
DO I=1,NFRZ
ITMP(IFRZ(I))=1
END DO
DO I=1,NBAS
IFRZ(I)=ITMP(I)
END DO
END IF
C
IF (NBAS.GT.NBASM) THEN
WRITE(6,*) 'NBAS, MAXIMUM: ',NBAS,NBASM
STOP 'INCREASE NBASM IN param.h AND RECOMPILE'
END IF
WRITE(6,*)
WRITE(6,*) ' SYSTEM: (and compiled maximum dimensions)'
WRITE(6,*) ' NATOMS ',NATOM,'(',NATMX,')'
WRITE(6,*) ' NBAS ',NBAS ,'(',NBASM,')'
WRITE(6,*) ' LMAX ',LMAX-1,' (',NLMAX-1,')'
WRITE(6,*)
WRITE(6,*)
WRITE(6,*) ' NUMBER OF SHELLS: ',NSHL
WRITE(6,*) ' NUMBER OF PRIMITIVES: ',IPRIM
WRITE(6,*) ' NUMBER OF BASIS FUNCTIONS:',NBAS
WRITE(6,*)
WRITE(6,*) ' calculating nuclear repulsion '
C
EN=0.D0
DO IAT=1,NATOM
XINUC=DBLE(NZ(IAT))
DO JAT=IAT+1,NATOM
RIJ=0.D0
DO I=1,3
RIJ=RIJ+(POS(I,IAT)-POS(I,JAT))*(POS(I,IAT)-POS(I,JAT))
END DO
EN=EN+XINUC*DBLE(NZ(JAT))/SQRT(RIJ)
END DO
END DO
C
WRITE(6,9901) EN
9901 FORMAT(//,' nuclear repulsion: ',F20.12,/)
C
C read overlap matrix
C
CALL RDMAT('OVERLAP',SAO)
WRITE(6,9902) (I,SAO(I,I),I=1,NBAS)
IF (L6D) THEN
DO I=1,NBAS
IF (ABS(SAO(I,I)).LE.1.D-2) THEN
WRITE(6,*) ' your overlap matrix seems to be'
$ ,' created with 5 d functions'
WRITE(6,*) ' please correct input to DALTON '
STOP ' 5D instead of 6D? '
END IF
END DO
END IF
9902 FORMAT(/,' THE DIAGONAL ELEMENTS OF THE OVERLAP MATRIX',/,50(4(I5
$ ,F10.3),/))
WRITE(6,*)
WRITE(6,*) ' READ OVERLAP MATRIX INTEGRALS '
C
C read Hamilton matrix
C
CALL RDMAT('HAMILTO',HAM)
WRITE(6,*) ' READ HAMILTON INTEGRALS '
C
C read additional one-electron integral matrix
C
IF (LMULTP) THEN
C we use the array for the kinetic energy integrals
CALL RDMAT('MULTINT',XMULP)
WRITE(6,*) ' READ MULTIPOLE INTEGRALS '
DO IALPH=1,NBAS
DO IBETA=1,NBAS
HAM(IALPH,IBETA)=HAM(IALPH,IBETA)+XMULP(IALPH,IBETA)
END DO
END DO
END IF
C
C read kinetic energy matrix
C
CALL RDMAT('KINETIC',XKIN)
WRITE(6,*) ' READ KINETIC ENERGY INTEGRALS '
WRITE(6,*)
C we have to substract the kinetic energy from the Hamilton integrals
DO I=1,NBAS
DO J=1,NBAS
HAM(I,J)=HAM(I,J)-XKIN(I,J)
END DO
END DO
C
C read startup vector
C
C in the canonical case we need just the occupation numbers, but we
C read the input correctly
C
DO I=1,NBAS
DO J=1,NBAS
CI(I,J)=0.D0
END DO
CI(I,I)=1.D0
END DO
C
OPEN(UNIT=IUNITO,STATUS='OLD',FORM='FORMATTED',FILE='GUESS',ERR
$ =6100)
READ(IUNITO,'(A)') BBBB
IF (BBBB.EQ.'E') THEN
WRITE(6,*)
WRITE(6,*) ' USING THE UNIT MATRIX AS START-UP VECTOR'
WRITE(6,*)
ELSE
REWIND(IUNITO)
WRITE(6,*)
WRITE(6,*) ' USING THE EXPLICITELY GIVEN START-UP VECTOR'
WRITE(6,*)
WRITE(6,*)
DO IFUNC=1,NBAS
READ(IUNITO,*,IOSTAT=KK) (CI(I,IFUNC),I=1,NBAS)
IF (KK.NE.0) THEN
WRITE(6,*)
WRITE(6,*) ' ERROR IN FUNCTION ',IFUNC
WRITE(6,*) ' LAST FUNCTION READ:'
WRITE(6,'(9F8.4)') (CI(I,IFUNC-1),I=1,NBAS)
STOP ' CORRECT INPUT!'
END IF
END DO
END IF
C
READ(IUNITO,*) (IOCC(I),I=1,NBAS)
READ(IUNITO,*) (IOCCS(I),I=1,NBAS)
CLOSE(IUNITO)
C
C for the moment, only closed shells
C
NOCC=0
DO I=1,NBAS
IF (IOCC(I).NE.0.AND.IOCC(I).NE.IOCCS(I))
- STOP 'NOT ALL SHELLS CLOSED'
IF (IOCC(I).EQ.IOCCS(I)) THEN
NOCC=NOCC+1
END IF
END DO
WRITE(6,*) ' NOCC ',NOCC
IF (LREMO) THEN
IF (NOCC.NE.NREMO) THEN
WRITE(6,*) ' your NREMO seems not correct '
LREMO=.FALSE.
END IF
END IF
WRITE(6,*)
NVIRT =NBAS-NOCC
C
IF (LTST) STOP ' YOU ASKED FOR A TEST RUN '
C
CALL TIMING('READ')
C
IF (.NOT.LCOREH) THEN
C sorting of the occupied to the left
WRITE(6,'(30I2)') (IOCC (I),I=1,NBAS)
WRITE(6,'(30I2)') (IOCCS(I),I=1,NBAS)
CALL VSRTO(CI,IOCC,IOCCS)
WRITE(6,'(30I2)') (IOCC (I),I=1,NBAS)
WRITE(6,'(30I2)') (IOCCS(I),I=1,NBAS)
CALL TIMING('VRST')
END IF
C
C now we start to work
C the SCF
C
IF (LCAN) THEN
CALL PRECAN
CALL SCFCAN
CALL TIMING('SCF ')
CALL WVEC(2)
ELSE
C
CALL SCFLOC
CALL TIMING('SCF ')
CALL WVEC(1)
END IF
C
C we are done with the SCF loops, we may calculate a Fock
C matrix corresponding to the orbitals, without further treatment
CALL CFOCKN
CALL TIMING('WAVE')
C
IF (LPIPM) THEN
CALL PIPEKM
CALL TIMING('P-M ')
CALL WVEC(3)
END IF
C
IF (LBOYS) THEN
CALL BOYS
CALL TIMING('BOYS')
CALL WVEC(4)
END IF
C
IF (LQMC) THEN
CALL OUTQMC
CALL TIMING('QMC ')
END IF
C
IF (LWFN) THEN
C CALL OUTWFN
WRITE(6,*) ' Gaussian WFN not yet implemented, work harder '
CALL TIMING('WFU ')
END IF
C
C calculate dipole moment
IF (LDIPOL) THEN
CALL DIPOLC
CALL TIMING('DIPL')
CALL QPOLC
CALL TIMING('QUAD')
END IF
C
IF (LMOLD) THEN
CALL MOLDEN
CALL TIMING('MOLD')
END IF
C
IF (LICSCF) THEN
CALL ICSCF
CALL TIMING('ICSF')
END IF
C
X=CPTIME(4)
CALL DATING(PNAME,2)
STOP
901 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
6100 CONTINUE
WRITE(6,*) ' NO FILE as starting vector FOUND ... '
END
C
SUBROUTINE SCFLOC
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
LOGICAL LCONV
C
THRESH=1.D-10
EOLD=0.D0
ITER=0
LCONV=.FALSE.
LMIX=.FALSE.
C
100 CONTINUE
CALL ORTHO(CI,NOCC,NBAS)
CALL TIMING ('ORTH')
C
C construct density matrix
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
SUM=0.D0
DO I=1,NOCC
SUM=SUM+CI(IALPH,I)*CI(IBETA,I)
END DO
P(IALPH,IBETA)=2.D0*SUM
END DO
END DO
C
C construct Fock matrix, and calculate the total energy
C
CALL CFOCK
CALL TIMING ('FOCK')
C
C
IF (ITER.GT.5) LMIX=.TRUE.
IF (LEXTER.AND.ITER.GT.2) LMIX=.TRUE.
IF (LMIX) THEN
C save the Fock matrix in atomic orbitals
OPEN(UNIT=15,FILE='FOCK.TMP',STATUS='UNKNOWN',FORM='FORMATTED')
DO I=1,NBAS
DO J=1,NBAS
WRITE(15,*) F(I,J)
END DO
END DO
CLOSE(15)
END IF
C
C transform the Fock matrix in molecular orbitals
CALL TRANS(F,CI,NBAS,NBAS)
C
C reorder the vectors according to the orbital energy
DO I=1,NBAS
ORBEN(I)=F(I,I)
END DO
c$$$C
c$$$C now we sort quite a lot at a time: ORBEN, CI, and F,
c$$$C by BUBBLESORT, despite the lecture of the Numerical Recipes
c$$$C
c$$$ CALL BUBBLE(NBAS,NBAS,ORBEN,F,CI)
c$$$C
WRITE(6,9901)
DO I=1,NBAS
WRITE(6,9902) I,ORBEN(I)
END DO
9901 FORMAT(' ORBITAL ENERGIES: ',/,' I ORBEN',/
$ ,' =========================================')
9902 FORMAT(I7,F14.6)
WRITE(6,*)
WRITE(6,*)
C
C sum the occ-virt matrix elements of the Fock matrix
SF=0.D0
DO I=1,NOCC
DO J=NOCC+1,NBAS
SF=SF+F(I,J)*F(I,J)
END DO
END DO
SF=SQRT(SF)/DBLE(NOCC*(NBAS-NOCC))
C
WRITE(6,9903) ITER,ETOT,ETOT-EOLD,SF
9903 FORMAT(' ITER',I4,' EA=',F12.4,' DEA=',E8.2,' SF=',E8.2)
C
IF (LECONV) THEN
C convergence with respect to energy
IF (ABS(ETOT-EOLD).LE.THRESH) LCONV=.TRUE.
ELSE
C convergence with respect to SFA and SFB
IF (SF.LE.THRESH) LCONV=.TRUE.
END IF
IF (LCONV) GO TO 200
C
IF (LCICOM) THEN
C we transform the bielectronic integrals for the complete SCI matrix
CALL TRANSDA
CALL TIMING('TRAN')
END IF
C
CALL DIAGFO(F,CI,NOCC,NBAS)
CALL TIMING ('DIAG')
C
EOLD=ETOT
ITER=ITER+1
IF (ITER.GE.NITMAX) GO TO 300
GO TO 100
C
200 CONTINUE
C we return
WRITE(6,9905) ITER
9905 FORMAT(' CONVERGENCE AFTER ',I4
$ ,' ITERATIONS ')
RETURN
C too many iterations
300 CONTINUE
WRITE(6,9906) NITER
9906 FORMAT(' WE USED ALL OF THE ',I4,' ITERATIONS; NO CONVERGENCE ')
RETURN
END
C
SUBROUTINE SCFCAN
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
LOGICAL LCONV
C
C just plain canonical SCF
C
IF (LCOREH) THEN
C arrange IOCC
DO I=1,NOCC
IOCC(I)=2
END DO
DO I=NOCC+1,NBAS
IOCC(I)=0
END DO
ITER=0
ELSE
ITER=1
END IF
C
EOLD=0.D0
LCONV=.FALSE.
C LECONV=.TRUE.
THRESH=1.D-10
LMIX=.FALSE.
C
100 CONTINUE
C
C construct density matrix
C
IF (ITER.NE.0) THEN
DO IALPH=1,NBAS
DO IBETA=IALPH,NBAS
SUM=0.D0
DO I=1,NOCC
SUM=SUM+CI(IALPH,I)*CI(IBETA,I)
END DO
P(IALPH,IBETA)=2.D0*SUM
END DO
END DO
C
DO IALPH=1,NBAS
DO IBETA=1,IALPH-1
P(IALPH,IBETA)=P(IBETA,IALPH)
END DO
END DO
C
C construct Fock matrix, and calculate the total energy
C
CALL CFOCK
ELSE
C for the first iteration we take just the core Hamiltonian
WRITE(6,*)
$ ' FOR THE FIRST ITERATION WE DIAGONALIZE'
$ ,' THE CORE HAMILTONIAN '
DO IALPH=1,NBAS
DO IBETA=IALPH+1,NBAS
F(IALPH,IBETA)=HAM(IALPH,IBETA)+XKIN(IALPH,IBETA)
END DO
C starting with zero on the diagonal seems a better start (?)
F(IALPH,IALPH)=0.D0
C this is the real core hamiltonian
C F(IALPH,IALPH)=HAM(IALPH,IALPH)+XKIN(IALPH,IALPH)
END DO
DO IALPH=1,NBAS
DO IBETA=IALPH+1,NBAS
F(IBETA,IALPH)=F(IALPH,IBETA)
END DO
END DO
C
END IF
C
IF (ITER.GT.5) LMIX=.TRUE.
IF (LEXTER.AND.ITER.GT.2) LMIX=.TRUE.
IF (LMIX) THEN
C save the Fock matrix in atomic orbitals
OPEN(UNIT=15,FILE='FOCK.TMP',STATUS='UNKNOWN',FORM='FORMATTED')
DO I=1,NBAS
DO J=1,NBAS
WRITE(15,*) F(I,J)
END DO
END DO
CLOSE(15)
END IF
C
C transform the Fock matrix in molecular orbitals
CALL TRANS(F,CI,NBAS,NBAS)
C
C reorder the vectors according to the orbital energy
DO I=1,NBAS
ORBEN(I)=F(I,I)
END DO
C
C now we sort quite a lot at a time: ORBEN, CI, and F,
C by BUBBLESORT, we do not know how to do more elegant
C
CALL BUBBLE(NBAS,NBAS,ORBEN,F,CI)
C
c$$$ WRITE(6,9901)
c$$$ DO I=1,NBAS
c$$$ WRITE(6,9902) I,ORBEN(I)
c$$$ END DO
c$$$ 9901 FORMAT(' ORBITAL ENERGIES: ',/
c$$$ $ ,' =========================')
c$$$ 9902 FORMAT(I7,F14.6)
c$$$ WRITE(6,*)
WRITE(6,*)
C
C sum the occ-virt matrix elements of the Fock matrix
SFN=0.D0
SUMO=0.D0
DO I=1,NOCC
SUMO=SUMO+ORBEN(I)
DO J=NOCC+1,NBAS
SFN=SFN+F(I,J)*F(I,J)
END DO
END DO
SFN=SQRT(SFN)/DBLE(NOCC*(NBAS-NOCC))
C
WRITE(6,9903) ITER,ETOT,ETOT-EOLD,SFN
9903 FORMAT(' ITER ',I4,': E_TOT=',F20.12,' EDIFF=',E8.2,' SFN=',E8.2)
WRITE(6,*) ' sum of orbital energies : ',2.D0*SUMO
c$$$C construct again the HF energy from HAM and ORBEN
c$$$ CALL RDMAT('HAMILTO',HAM)
c$$$ CALL TRANS(HAM,CI,NBAS,NBAS)
c$$$ EHF=EN
c$$$ DO I=1,NOCC
c$$$ EHF=EHF+HAM(I,I)+ORBEN(I)
c$$$ END DO
c$$$ WRITE(6,*) ' reconstructed HF energy ',EHF
c$$$
c$$$ CALL RDMAT('HAMILTO',HAM)
C
IF (ITER.NE.0) THEN
IF (LECONV) THEN
C convergence with respect to energy
IF (ABS(ETOT-EOLD).LE.THRESH)
$ LCONV=.TRUE.
ELSE
C convergence with respect to SFA and SFB
IF (SFN.LE.THRESH) LCONV=.TRUE.
END IF
IF (LCONV) GO TO 200
END IF
C
C diagonalize the Fock matrix
CALL GENCAN(F,CI,ORBEN,NBAS)
C we save the orbitals for a possible restart
CALL WVEC(5)
C
EOLD=ETOT
ITER=ITER+1
IF (ITER.GT.NITMAX) GO TO 300
GO TO 100
C
200 CONTINUE
C we return
WRITE(6,9905) ITER
9905 FORMAT(' CONVERGENCE AFTER ',I4,' ITERATIONS')
RETURN
C too many iterations
300 CONTINUE
WRITE(6,9906) NITMAX
9906 FORMAT(' WE USED ALL OF THE ',I4,' ITERATIONS; NO CONVERGENCE')
RETURN
END
C
SUBROUTINE ORTHO(CI,NOCC,NBAS)
C we do not include the common 'common_syst.h' here
C we do not include the common 'common_vect.h'
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION CI(NBASM,NBASM)
C
CALL SHALF(CI,NOCC,NBAS)
C
CALL SFULL(CI,NOCC,NBAS)
C
CALL SHALF(CI(1,NOCC+1),NBAS-NOCC,NBAS)
C
C CALL PFUNC(CI,NBAS,NBAS)
RETURN
END
C
SUBROUTINE SNORM(CI,NVECT,NBAS)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION CI(NBASM,NVECT)
C
C normalisation of vectors
C
DO I=1,NVECT
SL=0.D0
DO K=1,NBAS
CNN=0.D0
DO L=1,NBAS
CNN=CNN+CI(L,I)*SAO(K,L)
END DO
SL=SL+CI(K,I)*CNN
END DO
SL=1.D0/SQRT(SL)
DO J=1,NBAS
CI(J,I)=CI(J,I)*SL
END DO
C WRITE(6,*) ' VECTOR: ',I,' NORM=',SL
END DO
C
RETURN
END
C
SUBROUTINE SHALF(CI,NOCC,NBAS)
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION SMOTR(NBASM*(NBASM+1)/2)
DIMENSION CWRK(NBASM,NBASM),EW(NBASM),C3(NBASM),C4(NBASM)
DIMENSION CWRK2(NBASM,NBASM)
DIMENSION CI(NBASM,NBASM)
DO I=1,NBAS
DO J=1,NBAS
CWRK(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(CWRK,CI,NOCC,NBAS)
OPEN(UNIT=13,FILE='SMO.TMP',FORM='FORMATTED',STATUS='UNKNOWN')
DO I=1,NBAS
DO J=1,I
IF (ABS(CWRK(I,J)).GT.1.D-9) WRITE(13,'(2I5,F20.12)') I,J,CWRK(I
$ ,J)
END DO
END DO
CLOSE(13)
IF (.NOT.LCOREH) THEN
DO I=1,NOCC
IF (ABS(1.D0-CWRK(I,I)).GT.1.D-3) THEN
WRITE(6,9007) I,CWRK(I,I)
9007 FORMAT(' bad norm? DIAGONAL ELEMENT No ',I4,': ',F15.8)
CALL PFUNC(CI(1,I),1,NBAS)
END IF
END DO
END IF
INDX=1
DO I=1,NOCC
DO J=1,I
SMOTR(INDX)=CWRK(I,J)
INDX=INDX+1
END DO
END DO
C
C diagonalisation of SMOTR
C
NM=NBASM
NS=NOCC
NV=NOCC*(NOCC+1)/2
CALL RSP(NBASM,NOCC,NV,SMOTR,EW,1,CWRK,C3,C4,IERR)
WRITE(6,*)
WRITE(6,*) ' SHALF: SMALLEST EIGENVALUE',EW(1)
WRITE(6,*)
C DO I=1,NOCC
C WRITE(6,*) ' EIGENWERT ',I,EW(I)
C WRITE(6,'(4(I4,F12.6))') (J,CWRK(J,I),J=1,NBAS)
C END DO
C
IF (EW(1).LT.EPS) THEN
WRITE(6,*) ' WE HAVE EPS = ' ,EPS,' AND EW(1) = ',EW(1)
STOP 'SHALF: LIN DEPENDENCIES'
END IF
C
DO I=1,NOCC
DO J=1,NOCC
CWRK2(I,J)=0.D0
END DO
END DO
DO L=1,NOCC
SSS=1.D0/SQRT(EW(L))
DO J=1,NOCC
DO I=1,NOCC
CWRK2(I,J)=CWRK2(I,J)+SSS*CWRK(I,L)*CWRK(J,L)
END DO
END DO
END DO
C
C WRITE(6,'(2I4,F12.6)') ((I,J,CWRK2(I,J),I=1,NOCC),J=1,NOCC)
C
C the new functions
C
DO I=1,NOCC
DO K=1,NBAS
CWRK(K,I)=CI(K,I)
CI(K,I)=0.D0
END DO
END DO
DO I=1,NOCC
DO J=1,NOCC
CIJK=CWRK2(I,J)
DO K=1,NBAS
CI(K,I)=CI(K,I)+CIJK*CWRK(K,J)
END DO
END DO
END DO
C
C CALL PFUNC(CI,NOCC,NBAS)
RETURN
END
C
SUBROUTINE SFULL(CI,NOCC,NBAS)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
DIMENSION SMO(NBASM,NBASM),CIJK(NBASM)
DIMENSION CI(NBASM,NBASM)
NVIRT=NBAS-NOCC
DO I=1,NBAS
DO J=1,NBAS
SMO(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(SMO,CI,NBAS,NBAS)
C
C the corrected functions
C
ITER=0
NOV=NOCC+1
100 CONTINUE
DO I=NOV,NBAS
DO J=1,NOCC
CIJK(J)=SMO(I,J)
END DO
DO K=1,NBAS
CCC=0.D0
DO J=1,NOCC
CCC=CCC+CI(K,J)*CIJK(J)
END DO
CI(K,I)=CI(K,I)-CCC
END DO
END DO
C
DO I=1,NBAS
DO J=1,NBAS
SMO(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(SMO,CI,NBAS,NBAS)
C
C calculate new overlaps
SSS=0.D0
ITER=ITER+1
DO I=NOV,NBAS
DO J=1,NOCC
SSS=SSS+SMO(I,J)*SMO(I,J)
END DO
END DO
SSS=SQRT(SSS)/(NVIRT*NOCC)
WRITE(6,*) ' ITERATION, SSS:',ITER,SSS
IF (SSS.GT.EPS.AND.ITER.LE.100) GO TO 100
C
RETURN
END
C
SUBROUTINE CFOCK
INCLUDE 'param.h'
C
C construction of the Fock matrix with the ODA algorithm of E.Cancès included
C
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
C we have the "old" best density and the new density
C the density is PODA
COMMON /ODACOM/ PODA(NBASM,NBASM),XIODA(NBASM,NBASM),EMODA,EKODA
$ ,EHODA,EMPODA,LODAU
LOGICAL LODAU
C.. INCLUDE 'common_oda.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /UNITS/ IUNITR,IUNITX,IUNITO
C.. INCLUDE 'common_units.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
LOGICAL LREAD
DIMENSION XI(NBASM,NBASM)
C we construct the Fockmatrix, by reading only once
C the integral file
C
C optionally, we transform the (oo|vv) and (ov|ov) bielectronic
C integrals for the complete Singles-CI matrix
C
C count number of electrons again
c$$$ XEL=0.D0
c$$$ DO I=1,NBAS
c$$$ DO J=1,NBAS
c$$$ XEL=XEL+P(I,J)*SAO(I,J)
c$$$ END DO
c$$$ END DO
CALL DTRACE(SAO,P,NBAS,XEL)
WRITE(6,*)
WRITE(6,*) ' NUMBER OF ELECTRONS: expected:',NOCC*2
$ ,'; found: ',XEL
WRITE(6,*)
C
C the monoelectronic part
C
EONE=0.D0
EKIN=0.D0
EHAM=0.D0
ETWO=0.D0
DO IALPH=1,NBAS
DO IBETA=1,NBAS
XI(IALPH,IBETA)=0.D0
F(IALPH,IBETA)=XKIN(IALPH,IBETA)+HAM(IALPH,IBETA)
EONE=EONE+P(IALPH,IBETA)*F(IALPH,IBETA)
EKIN=EKIN+P(IALPH,IBETA)*XKIN(IALPH,IBETA)
EHAM=EHAM+P(IALPH,IBETA)*HAM(IALPH,IBETA)
END DO
END DO
IF (LMULTP) THEN
C
C the multipoles are already in the Hamilton integrals, so for calculating
C its contributions, we have to do this explicitely
C
EMULP=0.D0
DO IALPH=1,NBAS
DO IBETA=1,NBAS
EMULP=EMULP+P(IALPH,IBETA)*XMULP(IALPH,IBETA)
END DO
END DO
EHAM=EHAM-EMULP
END IF
C
C for the external construction we dump the vector onto VEC.TMP
C
IF (LEXTER) THEN
C
WRITE(6,*)
WRITE(6,*) ' dumping the vector in MOLPRO format onto VEC.TMP '
WRITE(6,*)
C
IUNITP=23
OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
- FILE='VEC.TMP')
WRITE(IUNITP,9234)
9234 FORMAT('# Orbitals created by ORTHO')
DO IALPH=1,NBAS
WRITE(IUNITP,'(3(E25.15,1H,))')
- (DBLE(ISIG(IPERM(IALPH)))*CI(IPERM(IALPH),J),J=1,NBAS)
END DO
CLOSE(IUNITP)
C
C call MOLPRO
C
WRITE(6,*)
WRITE(6,*) ' Calling Molpro ................. '
WRITE(6,*)
CALL SYSTEM
$ ('/home/reinh/bin/molpro02.7/molpro < MOLDFT.inp >> molpro.o
$ut')
C folded 1 (fixf $Revision: 1.3 $)
WRITE(6,*)
WRITE(6,*) ' Molpro terminated '
WRITE(6,*)
C
C read the Fock matrix
C
IUNITP=23
OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
- FILE='gorch.mat')
READ(IUNITP,*)
C we have to pay attention to the d functions, again!!!
READ(IUNITP,*) ((F(IPERM(IALPH),IPERM(IBETA))
- ,IALPH=1,NBAS),IBETA=1,NBAS)
CLOSE(IUNITP)
C and to the sign of f0, f3
DO IALPH=1,NBAS
XMULTA=ISIG(IALPH)
DO IBETA=1,NBAS
XMULTB=ISIG(IBETA)
F(IALPH,IBETA)=XMULTA*XMULTB*F(IALPH,IBETA)
END DO
END DO
OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
- FILE='dftfun.tmp')
READ(IUNITP,*)
READ(IUNITP,*)
C read the DFT energy
READ(IUNITP,*) COUENER,EXENER,DFTENER
CLOSE(IUNITP)
C
C we cannot calculate ETWO since we do not know the functional
ETWO=COUENER+XFAC*EXENER+DFTENER
C
ELSE
C
C the bielectronic part
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
XI(IALPH,IBETA)=0.D0
END DO
END DO
IUNITR=4
CALL WOPEN(IUNITR,2)
100 CONTINUE
CALL WREAD(LREAD)
IF (LAMBDA) THEN
DO III=1,INTRD
HBUF(III)=HBUF(III)*XLAMB
END DO
END IF
C//SKIPCI
cI XNU=2.D0
cI/ENDCI
DO III=1,INTRD
I=IBUF(III)
J=JBUF(III)
K=KBUF(III)
L=LBUF(III)
HHH=HBUF(III)
C
C//SKIPCI
cI CALL ACCUG(I,J,K,L,HHH,XNU,P,XI)
cI/ENDCI
C
XI(K,L)=XI(K,L)+4.D0*P(I,J)*HHH
XI(I,J)=XI(I,J)+4.D0*P(K,L)*HHH
XI(I,K)=XI(I,K)-P(J,L)*HHH
XI(J,K)=XI(J,K)-P(I,L)*HHH
XI(I,L)=XI(I,L)-P(J,K)*HHH
XI(J,L)=XI(J,L)-P(I,K)*HHH
C
END DO
IF (.NOT.LREAD) GO TO 100
C
CALL WCLOS(2)
C
C double the matrix ... for the non-diagonal elements
DO IALPH=1,NBAS
DO IBETA=IALPH+1,NBAS
XDUM=(XI(IALPH,IBETA)+XI(IBETA,IALPH))*0.5D0
XI(IALPH,IBETA)=XDUM
XI(IBETA,IALPH)=XDUM
END DO
END DO
C
IF (LODA) THEN
IF (LODAU) THEN
LODAU=.FALSE.
WRITE(6,*) ' first iteration, no new LAMBDA '
XLAMBD=1.D0
ELSE
C
C determine the optimal mixing parameter
C
A2=EONE
A1=EMODA
B1=0.D0
B2=0.D0
B3=0.D0
DO IALPH=1,NBAS
DO IBETA=1,NBAS
B3=B3+P(IALPH,IBETA) *XI(IALPH,IBETA)
B2=B2+P(IALPH,IBETA) *XIODA(IALPH,IBETA)
B1=B1+PODA(IALPH,IBETA)*XIODA(IALPH,IBETA)
END DO
END DO
B1=0.5D0*B1
B3=0.5D0*B3
B=A2-A1-2.D0*B1+B2
A=B1-B2+B3
WRITE(6,*) ' ODA: ONEOLD = ',A1
WRITE(6,*) ' ODA: ONENEW = ',A2
WRITE(6,*) ' ODA: PGOO = ',B1
WRITE(6,*) ' ODA: PGON = ',B2
WRITE(6,*) ' ODA: PGNN = ',B3
WRITE(6,*) ' ODA: A=',A,' B= ',B
IF (A.LT.0.D0) THEN
WRITE(6,*) ' negative quadratic coefficient !! '
XLAMBD=2.D0
ELSE
IF (ABS(A).GT.1.D-12) THEN
XLAMBD=-0.5D0*B/A
ELSE
XLAMBD=0.5D0
END IF
END IF
WRITE(6,*) ' best mixing parameter LAMBDA = ',XLAMBD
IF (XLAMBD.GT.1.D0) THEN
WRITE(6,*) ' setting LAMBDA = 1 '
XLAMBD=1.D0
END IF
ETOTOL=A1+B1+EN
ETOTNE=A2+B3+EN
ETOTMP=EN+A1+B1+XLAMBD*(B+XLAMBD*A)
WRITE(6,*) ' old energy : '
WRITE(6,*) ' total : ',ETOTOL
WRITE(6,*) ' mono : ',A1
WRITE(6,*) ' bi : ',B1
WRITE(6,*) ' new energy : '
WRITE(6,*) ' total : ',ETOTNE
WRITE(6,*) ' mono : ',A2
WRITE(6,*) ' bi : ',B3
WRITE(6,*) ' best new energy : ',ETOTMP
END IF
C
C construct the new quantities
ONELAM=1.D0-XLAMBD
C new best monoelectronic energies
C total mono
EMODA =ONELAM*EMODA+XLAMBD*EONE
C kinetic
EKODA =ONELAM*EKODA+XLAMBD*EKIN
C Hamilton (without multipoles)
EHODA =ONELAM*EHODA+XLAMBD*EHAM
C multipoles
EMPODA=ONELAM*EMPODA+XLAMBD*EMULP
C new PODA and XIODA matrix
DO IALPH=1,NBAS
DO IBETA=1,NBAS
PODA(IALPH,IBETA)=ONELAM*PODA(IALPH,IBETA)
- +XLAMBD*P(IALPH,IBETA)
XIODA(IALPH,IBETA)=ONELAM*XIODA(IALPH,IBETA)
- +XLAMBD*XI(IALPH,IBETA)
END DO
END DO
C new Fock matrix : F contains only the monoelectronic part of
C the integrals, no contraction of integrals with a density matrix
C
C we can add simply the new XIODA matrix to F
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
F(IALPH,IBETA)=F(IALPH,IBETA)
- +XIODA(IALPH,IBETA)
ETWO=ETWO+PODA(IALPH,IBETA)*XIODA(IALPH,IBETA)
END DO
END DO
EONE=EMODA
EKIN=EKODA
EHAM=EHODA
EMULP=EMPODA
C
ELSE
C add the simple bielectronic part to the Fock matrix
DO IALPH=1,NBAS
DO IBETA=1,NBAS
F(IALPH,IBETA)=F(IALPH,IBETA)+XI(IALPH,IBETA)
ETWO=ETWO+P(IALPH,IBETA)*XI(IALPH,IBETA)
END DO
END DO
END IF
ETWO=0.5D0*ETWO
C
c$$$C write the Fock matrix in MOLPRO format onto file
c$$$C
c$$$ IUNITP=23
c$$$ OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
c$$$ - FILE='FOCKM.TMP')
c$$$ WRITE(IUNITP,9235)
c$$$ 9235 FORMAT('# Fock matrix constructed by ORTHO')
c$$$ DO IALPH=1,NBAS
c$$$ WRITE(IUNITP,'(5(F15.8,1H,))') (F(IPERM(IALPH),IPERM(IBETA))
c$$$ $ ,IBETA=1,NBAS)
c$$$ END DO
c$$$ CLOSE(IUNITP)
c$$$
END IF
C Fock matrix mixing
IF (LMIX) THEN
OPEN(FILE='FOCK.TMP',FORM='FORMATTED',STATUS='UNKNOWN',UNIT=15
$ ,ERR=123)
WRITE(6,*) ' FOCK MATRIX MIXING WITH ',XFMIX*100.0D0,' PERCENT'
DO I=1,NBAS
DO J=1,NBAS
READ(15,*) FFA
c$$$ WRITE(6,*) ' I,J, OLD, NEW ',I,J,FFA,FFB,FA(I,J),FB(I,J)
F(I,J)=(1.D0-XFMIX)*F(I,J)+XFMIX*FFA
END DO
END DO
CLOSE(15,STATUS='DELETE')
123 CONTINUE
END IF
C
C here we arrive after the construction of the Fock matrix
C
ETOT=EONE+ETWO+EN
C
C save the Fock matrix for the program vind
OPEN(FILE='FOCK',FORM='FORMATTED',STATUS='UNKNOWN',UNIT=15)
WRITE(15,*) NBAS
WRITE(15,'(4E20.11)') ((F(IALPH,IBETA),IALPH=1,NBAS),IBETA=1
$ ,NBAS)
WRITE(15,*)
WRITE(15,'(4E20.11)') ETOT,EN,EONE,ETWO
CLOSE(15)
C
C the total energy
C
IF (LMULTP) THEN
WRITE(6,9902) EKIN,EHAM,EMULP,EKIN+EHAM+EMULP,ETWO,ETOT-EN
9902 FORMAT(' THE TOTAL ENERGY: ',/
$ ,' kinetic ',F20.12,/
$ ,' el-nucl ',F20.12,/
$ ,' el-multipoles ',F20.12,/
$ ,' (one-el ',F20.12,')',/
$ ,' el-el ',F20.12,/
$ ,' =====================================================',/
$ ,' total ',F20.12,' (without nuclear energy!)' ,/)
ELSE
WRITE(6,9901) EN,EKIN,EHAM,EKIN+EHAM,ETWO,ETOT,ETOT/EKIN
9901 FORMAT(' THE TOTAL ENERGY: ',/
$ ,' nuclear ',F20.12,/
$ ,' kinetic ',F20.12,/
$ ,' el-nucl ',F20.12,/
$ ,' (one-el ',F20.12,')',/
$ ,' el-el ',F20.12,/
$ ,' =====================================================',/
$ ,' total ',F20.12,' (virial theorem: ',F11.8,')',/)
END IF
RETURN
END
C
SUBROUTINE LEXSRT(IH0,H0,N)
IMPLICIT NONE
INTEGER N,L,IR,I,J,III
INTEGER IH0(4,N)
INTEGER IRA(4)
DOUBLE PRECISION H0(N),H0R
LOGICAL LCOMP4,LACT
C
C SORT BY IH0
C
IF (N.LT.2) GO TO 200
L=N/2+1
IR=N
10 CONTINUE
IF(L.GT.1)THEN
L=L-1
DO III=1,4
IRA(III)=IH0(III,L)
END DO
H0R=H0(L)
ELSE
DO III=1,4
IRA(III)=IH0(III,IR)
IH0(III,IR)=IH0(III,1)
END DO
H0R=H0(IR)
H0(IR)=H0(1)
IR=IR-1
IF (IR.EQ.1) THEN
DO III=1,4
IH0(III,1)=IRA(III)
END DO
H0(1)=H0R
GO TO 200
END IF
END IF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF (J.LT.IR)THEN
LACT=LCOMP4(IH0(1,J),IH0(1,J+1))
IF (LACT) THEN
J=J+1
END IF
END IF
LACT=LCOMP4(IRA,IH0(1,J))
IF(LACT)THEN
DO III=1,4
IH0(III,I)=IH0(III,J)
END DO
H0(I)=H0(J)
I=J
J=J+J
ELSE
J=IR+1
END IF
GOTO 20
END IF
DO III=1,4
IH0(III,I)=IRA(III)
END DO
H0(I)=H0R
GOTO 10
C
200 CONTINUE
RETURN
END
C
FUNCTION LCOMP4(IARR,JARR)
LOGICAL LCOMP4
INTEGER IARR(4),JARR(4)
C
IF (IARR(1).LT.JARR(1)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(1).GT.JARR(1)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(2).LT.JARR(2)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(2).GT.JARR(2)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(3).LT.JARR(3)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(3).GT.JARR(3)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(4).LT.JARR(4)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(4).GT.JARR(4)) THEN
LCOMP4=.FALSE.
ELSE
STOP ' EQUAL INDICES ENCOUNTERED! '
END IF
END IF
END IF
END IF
RETURN
END
C
FUNCTION HFIND(I1,J1,K1,L1,H0,IH0,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION IH0(4,N),H0(N)
LOGICAL LUP,LI,LJ,LK,LL
C
I=I1
J=J1
K=K1
L=L1
IF (I.GT.J) THEN
IDUM=I
I=J
J=IDUM
END IF
IF (K.GT.L) THEN
IDUM=K
K=L
L=IDUM
END IF
IF (I.GT.K) THEN
IDUM=I
I=K
K=IDUM
IDUM=J
J=L
L=IDUM
END IF
IF (I.EQ.K.AND.J.GT.L) THEN
IDUM=J
J=L
L=IDUM
END IF
C
C WRITE(6,*) ' HFIND: ',I1,J1,K1,L1,I,J,K,L
C
IND=N/2
INC=(N+1)/4+2
C
100 CONTINUE
IF (IH0(1,IND).LT.I) THEN
LUP=.TRUE.
ELSE IF (IH0(1,IND).GT.I) THEN
LUP=.FALSE.
ELSE IF (IH0(1,IND).EQ.I) THEN
IF (IH0(2,IND).LT.J) THEN
LUP=.TRUE.
ELSE IF (IH0(2,IND).GT.J) THEN
LUP=.FALSE.
ELSE IF (IH0(2,IND).EQ.J) THEN
IF (IH0(3,IND).LT.K) THEN
LUP=.TRUE.
ELSE IF (IH0(3,IND).GT.K) THEN
LUP=.FALSE.
ELSE IF (IH0(3,IND).EQ.K) THEN
IF (IH0(4,IND).LT.L) THEN
LUP=.TRUE.
ELSE IF (IH0(4,IND).GT.L) THEN
LUP=.FALSE.
ELSE IF (IH0(4,IND).EQ.L) THEN
HFIND=H0(IND)
RETURN
END IF
END IF
END IF
END IF
C
c$$$ WRITE(6,*) ' HFIND: ',INC,LUP,':',(IH0(JJJ,IND),JJJ=1,4)
c$$$ $ ,' (',I,J,K,L,')'
C
IF (LUP) THEN
IND=MIN(IND+INC,N)
ELSE
IND=MAX(IND-INC,1)
END IF
INC=(INC+1)/2
C
c$$$ WRITE(6,*) ' HFIND, II: ',IND
C
IF (INC.EQ.1) THEN
LI=IH0(1,IND).EQ.I
LJ=IH0(2,IND).EQ.J
LK=IH0(3,IND).EQ.K
LL=IH0(4,IND).EQ.L
IF (LI.AND.LJ.AND.LK.AND.LL) THEN
HFIND=H0(IND)
RETURN
END IF
C +1?
IND=IND+1
LI=IH0(1,IND).EQ.I
LJ=IH0(2,IND).EQ.J
LK=IH0(3,IND).EQ.K
LL=IH0(4,IND).EQ.L
IF (LI.AND.LJ.AND.LK.AND.LL) THEN
HFIND=H0(IND)
RETURN
END IF
C -1?
IND = IND-2
LI=IH0(1,IND).EQ.I
LJ=IH0(2,IND).EQ.J
LK=IH0(3,IND).EQ.K
LL=IH0(4,IND).EQ.L
IF (LI.AND.LJ.AND.LK.AND.LL) THEN
HFIND=H0(IND)
RETURN
END IF
HFIND=0.D0
RETURN
END IF
GO TO 100
END
C
SUBROUTINE TRANS(A,CI,NOCC,NBAS)
INCLUDE 'param.h'
DIMENSION AJKLM(NBASM,NBASM)
DIMENSION A(NBASM,NBASM),CI(NBASM,NBASM)
C
DO J=1,NBAS
DO I=1,NBAS
AJKLM(I,J)=0.D0
END DO
END DO
C
DO J=1,NOCC
DO K=1,NBAS
BB=0.D0
DO L=1,NBAS
BB=BB+CI(L,J)*A(K,L)
END DO
AJKLM(K,J)=AJKLM(K,J)+BB
END DO
END DO
C
DO J=1,NBAS
DO I=1,NBAS
A(I,J)=0.D0
END DO
END DO
C
DO I=1,NOCC
C only J=I..NBAS, copying is cheaper
DO J=I,NOCC
SSS=0.D0
DO K=1,NBAS
SSS=SSS+AJKLM(K,J)*CI(K,I)
END DO
A(I,J)=A(I,J)+SSS
END DO
END DO
C
C copy A(I,J)=A(J,I)
C
DO I=1,NOCC
DO J=I+1,NOCC
A(J,I)=A(I,J)
END DO
END DO
C
RETURN
END
C
SUBROUTINE BUBBLE(NOCC,NBAS,ORBEN,F,CI)
INCLUDE 'param.h'
DIMENSION ORBEN(NOCC),CI(NBASM,NOCC),F(NBASM,NOCC)
C
C vectors are NBAS X NOCC
C F is NOCC X NOCC
C
DO I=1,NOCC-1
DO J=1,NOCC-I
IF (ORBEN(J).GT.ORBEN(J+1)) THEN
DO II=1,NBAS
CDUM=CI(II,J)
CI(II,J)=CI(II,J+1)
CI(II,J+1)=CDUM
END DO
DO II=1,NOCC
FDUM=F(II,J)
F(II,J)=F(II,J+1)
F(II,J+1)=FDUM
END DO
DO II=1,NOCC
FDUM=F(J,II)
F(J,II)=F(J+1,II)
F(J+1,II)=FDUM
END DO
ODUM=ORBEN(J)
ORBEN(J)=ORBEN(J+1)
ORBEN(J+1)=ODUM
END IF
END DO
END DO
RETURN
END
C
C encore the Davidson, in a specialized version to pass F to ATIMES
C
SUBROUTINE DAVID(F,NOCC,NBAS,N,AB,B,ENERGY,THR,MAXDAV,ITER,ERR)
INCLUDE 'param.h'
C
C DIAGONALIZATION SUBROUTINE FOR LARGE SYMMETRIC MATRICES
C E. R. DAVIDSON, J. COMP. PHYS 17, 87 (1975)
C
C ce morceau est du CASDI de Daniel MAYNAU
C
PARAMETER (NITM=50)
DIMENSION B(*),AB(*)
C the correction vector, could be put as pieces onto disk
DIMENSION Q(1+NBASM*NBASM/4),DIAG(1+NBASM*NBASM/4)
DIMENSION F(NBASM,NBASM)
C
DIMENSION BAB(NITM,NITM)
DIMENSION ALPHA(NITM,NITM),ALAM(NITM),P1(NITM),P2(NITM)
LOGICAL YW
WRITE(6,*) ' DAVIDSON: ',N,' X ',N,' MATRIX '
5 FORMAT (//10X,'EIGENVALUE NUMBER',I3)
10 FORMAT (3X,'ITER=',I3,' EIGENVALUE=',F17.12,
* ', CONVERGENCE =',D10.2)
11 FORMAT (1X,'DAV: ITER=',I3,' EIGENV=',E10.4,
* ', CONV =',D10.2)
20 FORMAT (///5X,'THE DAVIDSON PROCEDURE DOES NOT CONVERGE AFTER',
* I4,' ITERATIONS'///)
CF CALL FLUSH(6)
C
C fill the diagonal
CALL DICALC(F,NOCC,NBAS,DIAG)
C
C dump the matrix
c$$$ OPEN(UNIT=51,FILE='MATRIX',STATUS='UNKNOWN',FORM='FORMATTED')
c$$$ DO I=1,N
c$$$ WRITE(51,'(2I4,F20.12)') I,I,DIAG(I)
c$$$ DO J=1,N
c$$$ B(J)=0.D0
c$$$ END DO
c$$$ B(I)=1.D0
c$$$ CALL ATIMES(F,NOCC,NBAS,B,AB)
c$$$ DO J=1,I-1
c$$$ IF (ABS(AB(J)).GT.1.D-7) WRITE(51,'(2I4,F20.12)') I,J,AB(J)
c$$$ END DO
c$$$ END DO
c$$$ CLOSE(51)
c$$$ STOP ' DUMPED MATRIX'
C
ZERO=0.D0
ONE=1.D0
EPS=1.D-25
C
C we open 2 direct-access files
C
IUNIT1=76
OPEN(UNIT=IUNIT1,FILE='VECTOR.TMP',STATUS='UNKNOWN',ACCESS='DIRECT
$'
C folded 1 (fixf $Revision: 1.3 $)
$ ,RECL=N*8)
IUNIT2=87
OPEN(UNIT=IUNIT2,FILE='HVECT.TMP',STATUS='UNKNOWN',ACCESS='DIRECT'
$ ,RECL=N*8)
C
ITER=1
C
SUMQ=0.D0
DO I=1,N
SUMQ=SUMQ+B(I)*B(I)
END DO
SUMQ=SQRT(SUMQ)
SUMQQ=0.D0
DO I=1,N
AB(I)=0.D0
B(I)=B(I)/SUMQ
SUMQQ=SUMQQ+B(I)*B(I)
END DO
C
CALL ATIMES(F,NOCC,NBAS,B,AB)
C
CALL PUTV(B,1,N,IUNIT1)
CALL PUTV(AB,1,N,IUNIT2)
BABL=0.D0
DO J=1,N
BABL=BABL+B(J)*AB(J)
END DO
BAB(1,1)=BABL
C
M=1
C
C DIAGONALIZATION OF BAB ... EIGENVALUES ALAM, EIGENVECTORS ALPHA
C
C this is the iteration loop
C
70 CONTINUE
C
C this small matrix BAB has to be diagonalized
C
IERR=0
MATZ=1
CALL RS(NITM,M,BAB,ALAM,MATZ,ALPHA,P1,P2,IERR)
IF (IERR.NE.0) THEN
WRITE(6,*) ' ERROR IN DIAGONALIZATION OF Heff'
DO I=1,N
B(I)=0.D0
END DO
RETURN
END IF
C
C SETS UP THE CORRECTION VECTOR Q
C AND PERFORMS THE CONVERGENCE TEST
C
75 CONTINUE
C we need only the lowest eigenvalue
CALL GETV(B,1,N,IUNIT1)
CALL GETV(AB,1,N,IUNIT2)
AVAL=ALAM(1)
AVEC=ALPHA(1,1)
DO I=1,N
DIFFI=AB(I)-AVAL*B(I)
Q(I)=AVEC*DIFFI
END DO
DO L=2,M
CALL GETV(B,L,N,IUNIT1)
CALL GETV(AB,L,N,IUNIT2)
C ... and the corresponding eigenvector
AVEC=ALPHA(L,1)
DO I=1,N
DIFFI=AB(I)-AVAL*B(I)
Q(I)=Q(I)+AVEC*DIFFI
END DO
END DO
C
C convergence?
SUMQ=0.D0
DO I=1,N
SUMQ=SUMQ+Q(I)*Q(I)
END DO
SUMQ=DSQRT(SUMQ/N)
c$$$ WRITE (6,10) ITER,ALAM(1),SUMQ
c$$$ IF (ITER.EQ.1) CALL TIMING('IT 1')
c$$$ IF (ITER.EQ.2) CALL TIMING('IT 2')
CF CALL FLUSH(6)
IF (ITER.GT.2)THEN
IF (SUMQ.LT.THR) GO TO 200
END IF
IF (ITER.GE.MAXDAV) GO TO 250
ITER=ITER+1
C
C IF M OVERFLOWS, THE B SUBSPACE IS THOROUGHLY REDEFINED
C
IF (M.GE.NITM) GO TO 200
C
C COMPUTES THE VECTOR B(M+1) TO ENLARGE THE B SUBSPACE
C
C Q contains (H-E_0)B
C
IQ=0
DO 100 I=1,N
IQ=IQ+1
IF (DABS(Q(IQ))-EPS) 95,95,97
95 CONTINUE
Q(IQ)=ZERO
GO TO 100
97 CONTINUE
IF(DABS(ALAM(1)-DIAG(I)).LT.1.D-15) THEN
WRITE(6,*)'QUASI DEGENERESCENCE',ALAM(1),DIAG(I)
Q(IQ)=Q(IQ)/(ALAM(1)-DIAG(I)-0.01)
ELSE
Q(IQ)=Q(IQ)/(ALAM(1)-DIAG(I))
END IF
100 CONTINUE
C
C -- ORTHOGONALISATION DE B(M+1) AUX B(1)...B(M)
C
DO L=1,M
CALL GETV(B,L,N,IUNIT1)
SUM=ZERO
DO I=1,N
SUM=SUM+B(I)*Q(I)
END DO
DO I=1,N
Q(I)=Q(I)-SUM*B(I)
END DO
END DO
SUM=ZERO
DO I=1,N
SUM=SUM+Q(I)*Q(I)
END DO
SUM=ONE/DSQRT(SUM)
DO I=1,N
B(I)=SUM*Q(I)
END DO
C
M=M+1
CALL PUTV(B,M,N,IUNIT1)
DO I=1,N
AB(I)=0.D0
END DO
CALL ATIMES(F,NOCC,NBAS,B,AB)
C
C normal CI
C
CALL PUTV(AB,M,N,IUNIT2)
DO L=1,M
SUM=ZERO
CALL GETV(B,L,N,IUNIT1)
DO J=1,N
SUM=SUM+B(J)*AB(J)
END DO
BAB(M,L)=SUM
END DO
C
C BACK TO THE DIAGONALIZATION OF BAB
C
GO TO 70
C
C THE EIGENVECTORS OF BAB ARE BACK-TRANSFORMED
C IN THE ORIGINAL BASIS
C
C the eigenvector is constructed from the stored intermediates,
C and iterations are restarted, if not all iterations have been used
C
200 CONTINUE
C
C construct vector
C
DO I=1,N
B(I)=ZERO
END DO
DO LL=1,M
ALP=ALPHA(LL,1)
CALL GETV(AB,LL,N,IUNIT1)
DO I=1,N
B(I)=B(I)+AB(I)*ALP
END DO
END DO
C
C on convergence we may return
C
IF (SUMQ.LT.THR) THEN
CLOSE(IUNIT1,STATUS='DELETE')
CLOSE(IUNIT2,STATUS='DELETE')
ENERGY=ALAM(1)
RETURN
END IF
C
C otherwise we restart the iterations from the
C improved starting vector
C
CALL PUTV(B,1,N,IUNIT1)
C
C construct H*vector
C
DO I=1,N
B(I)=ZERO
END DO
DO LL=1,M
ALP=ALPHA(LL,1)
CALL GETV(AB,LL,N,IUNIT2)
DO I=1,N
B(I)=B(I)+AB(I)*ALP
END DO
END DO
CALL PUTV(B,1,N,IUNIT2)
C
WRITE (6,9221)
9221 FORMAT(' DAVIDSON ITERATIONS ARE RESTARTED WITH NEW VECTOR')
C
CALL GETV( B,1,N,IUNIT1)
CALL GETV(AB,1,N,IUNIT2)
BABL=0.D0
DO I=1,N
BABL=BABL+B(I)*AB(I)
END DO
C
M=1
BAB(1,1)=BABL
ALPHA(1,1)=ONE
ALAM(1)=BABL
GO TO 75
250 CONTINUE
WRITE (6,20) ITER
C
C the procedure did not converge, however we return the best vector
C found
C
DO I=1,N
B(I)=ZERO
END DO
DO LL=1,M
ALP=ALPHA(LL,1)
CALL GETV(AB,LL,N,IUNIT1)
DO I=1,N
B(I)=B(I)+AB(I)*ALP
END DO
END DO
ENERGY=ALAM(1)
CLOSE(IUNIT1,STATUS='DELETE')
CLOSE(IUNIT2,STATUS='DELETE')
C
RETURN
END
C
SUBROUTINE GETV(X,IVEC,NDIM,IUNIT)
INCLUDE 'param.h'
DIMENSION X(NDIM)
C WRITE(6,*) ' WE TRY TO READ RECORD No',IVEC,' FROM FILE ',IUNIT
READ(IUNIT,REC=IVEC) X
RETURN
END
C
SUBROUTINE PUTV(X,IVEC,NDIM,IUNIT)
INCLUDE 'param.h'
DIMENSION X(NDIM)
C WRITE(6,*) ' WE TRY TO WRITE RECORD No',IVEC,' TO FILE ',IUNIT
WRITE(IUNIT,REC=IVEC) X
RETURN
END
C
SUBROUTINE ATIMES(F,NOCC,NBAS,VECT1,VECT2)
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
DIMENSION VECT1(*),VECT2(*),F(NBASM,NBASM)
LOGICAL LREAD
C
C we calculate the effect of the Singles -CI matrix on the vector VECT1
C
C the Singles-CI matrix:
C order is (i,a)
C
C 0 sqrt(2)*F_ia .......
C
S2=SQRT(2.D0)
NVIRT=NBAS-NOCC
C
C first line, the diagonal and initialization of VECT2
VECT2(1)=0.D0
INDX=1
DO I=1,NOCC
DO IA=NOCC+1,NBAS
INDX=INDX+1
INDX2=KPOS(I,IA,NOCC,NBAS)
IF (INDX.NE.INDX2) WRITE(6,*) ' KPOS: ',I,IA,INDX,INDX2
FIA=F(I,IA)
FACT=FIA*S2
VECT2(INDX)=FACT*VECT1(1)
VECT2(1) =FACT*VECT1(INDX) + VECT2(1)
C
C the diagonal
C
FACT=F(IA,IA)-F(I,I)
VECT2(INDX)=VECT2(INDX)+FACT*VECT1(INDX)
C
END DO
END DO
C
C the rest
C
DO I=1,NOCC
DO IA=NOCC+1,NBAS
C INDX1 this is the position of (I,IA) in the vector
INDX1=KPOS(I,IA,NOCC,NBAS)
C
DO IB=NOCC+1,IA-1
C = F(a,b)
FACT=F(IA,IB)
C INDX2 this is the position of (I,IB) in the vector
INDX2=KPOS(I,IB,NOCC,NBAS)
VECT2(INDX1)=VECT2(INDX1)+VECT1(INDX2)*FACT
VECT2(INDX2)=VECT2(INDX2)+VECT1(INDX1)*FACT
END DO
C
C = -F(i,j)
C
DO J=1,I-1
FACT=-F(I,J)
C INDX2 this is the position of (J,IA) in the vector
INDX2=KPOS(J,IA,NOCC,NBAS)
VECT2(INDX1)=VECT2(INDX1)+VECT1(INDX2)*FACT
VECT2(INDX2)=VECT2(INDX2)+VECT1(INDX1)*FACT
END DO
END DO
END DO
C
IF (LCICOM) THEN
C
C general element = 2(ia|jb)-(ij|ab) with ia > jb
C
C we read elements from a file i,j,k,l, h
C
C the integrals are on unit 12, we just read them
IUNITR=12
CALL WOPEN(IUNITR,2)
100 CONTINUE
CALL WREAD(LREAD)
DO III=1,INTRD
I =IBUF(III)
IA=JBUF(III)
J =KBUF(III)
IB=LBUF(III)
HHH=HBUF(III)
C
INDX1=KPOS(I,IA,NOCC,NBAS)
INDX2=KPOS(J,IB,NOCC,NBAS)
FACT=HHH
C INDX2 this is the position of (I,IB) in the vector
VECT2(INDX1)=VECT2(INDX1)+VECT1(INDX2)*FACT
VECT2(INDX2)=VECT2(INDX2)+VECT1(INDX1)*FACT
END DO
IF (.NOT.LREAD) GO TO 100
C we may even delete the file
CALL WCLOS(2)
END IF
RETURN
END
FUNCTION KPOS(I,IA,NOCC,NBAS)
IMPLICIT NONE
INTEGER I,IA,KPOS,NOCC,NBAS
KPOS=1+(NBAS-NOCC)*I+IA-NBAS
RETURN
END
C
SUBROUTINE DENS(CI)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /UNITS/ IUNITR,IUNITX,IUNITO
C.. INCLUDE 'common_units.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
DIMENSION CI(NBASM,NBASM)
C
C calculate number of electrons
C
WRITE(6,*) 'BASIS FUNCTION ASSIGNED NUMBER OF ELECTRONS'
FTOT=0.D0
DO K=1,NBAS
FF=0.D0
DO L=1,NBAS
FF=FF+P(K,L)*SAO(K,L)
END DO
FF=FF
FTOT=FTOT+FF
IF (ABS(FF).GT.1.D-6) WRITE(6,'(5X,I6,15X,F13.6)') K,FF
END DO
WRITE(6,*)
WRITE(6,*) 'TOTAL NUMBER OF ELECTRONS',FTOT
WRITE(6,*)
RETURN
END
C
SUBROUTINE DIAGFO(F,CI,NOCC,NBAS)
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
DIMENSION F(NBASM,NBASM),CI(NBASM,NBASM)
DIMENSION CIN(NBASM)
DIMENSION VECT1(1+NBASM*NBASM/4),VECT2(1+NBASM*NBASM/4)
C
C here is the cruxial step of the calculation: the diagonalization of
C the CI matrix; we use the Davidson procedure for this
C
NVIRT=NBAS-NOCC
NDIM=1+NOCC*NVIRT
THR=1.D-9
MAXDAV=250
C
VECT1(1)=1.D0
DO I=2,NDIM
VECT1(I)=0.D0
END DO
C
CALL DAVID(F,NOCC,NBAS,NDIM,VECT2,VECT1,ENERGY,THR,MAXDAV,NITERD,E
$RR)
C folded 1 (fixf $Revision: 1.3 $)
C
C VECT1 is the Eigenvector, ENERGY the lowering of the total energy
WRITE(6,*) ' CONVERGENCE after ',NITERD,' ITERATIONS'
WRITE(6,*) ' LOWERING of the total energy: ',ENERGY
WRITE(6,*) ' WEIGHT of the reference : ',VECT1(1)
C
C we have to construct new vectors
C
C occ: C(alpha,i) = C(alpha,i) + sum_a c_ia/c_1 * C(alpha,a)
C virt: C(alpha,a) = C(alpha,a) - sum_i c_ia/c_1 * C(alpha,i)
C
CNRM=VECT1(1)
WRITE(6,*) ' THE COEFFICIENTS OF THE EIGENVECTOR: '
WRITE(6,*) ' occ virt coeff '
WRITE(6,*) ' ============================= '
INDX=1
DO I=1,NOCC
DO IA=NOCC+1,NBAS
INDX=INDX+1
VECT1(INDX)=VECT1(INDX)/CNRM
IF (ABS(VECT1(INDX)).GT.1.D-4) WRITE(6,'(2I6,F20.12)') I,IA
$ ,VECT1(INDX)
END DO
END DO
WRITE(6,*)
C
INDX=1
DO I=1,NOCC
DO IA=NOCC+1,NBAS
INDX=INDX+1
IF (ABS(VECT1(INDX)).GT.0.5D0) THEN
VECDUM=VECT1(INDX)
VECT1(INDX)=SIGN(0.5D0,VECT1(INDX))
WRITE(6,*) ' SETTING COEFF ',I,' -> ',IA,' FROM ',VECDUM,' TO '
$ ,VECT1(INDX)
END IF
END DO
END DO
WRITE(6,*)
C
IF (LODA) THEN
XMIXX=1.D0
WRITE(6,*) ' ODA : coefficient mixing with 100% '
ELSE
XMIXX=0.6D0
WRITE(6,*) ' coefficient mixing with 60% '
END IF
C
DO IALPH=1,NBAS
C
DO I=1,NBAS
CIN(I)=0.D0
END DO
INDX=1
DO I=1,NOCC
DO IA=NOCC+1,NBAS
INDX=INDX+1
FACT=XMIXX*VECT1(INDX)
CIN(I) =CIN(I) + FACT*CI(IALPH,IA)
CIN(IA)=CIN(IA) - FACT*CI(IALPH,I)
END DO
END DO
C
DO I=1,NBAS
CI(IALPH,I)=CI(IALPH,I)+CIN(I)
END DO
C close loop over ialph
END DO
C
RETURN
END
C
SUBROUTINE DICALC(F,NOCC,NBAS,DIAG)
INCLUDE 'param.h'
DIMENSION F(NBASM,NBASM),DIAG(*)
C
C the diagonal of the Singles-CI matrix
C
c$$$ WRITE(6,*) ' DICALC: ',NOCC,NBAS
c$$$ WRITE(6,'(10H ORBITAL ,I4,F20.12)') (J,F(J,J),J=1,NBAS)
DIAG(1)=0.D0
INDX=1
DO I=1,NOCC
FII=F(I,I)
DO IA=1+NOCC,NBAS
INDX=INDX+1
DIAG(INDX)=F(IA,IA)-FII
END DO
END DO
RETURN
END
SUBROUTINE PFUNC(CI,NVECT,NBAS)
INCLUDE 'param.h'
DIMENSION CI(NBASM,NBASM),CIOUT(5)
C
DO IVEC=1,NVECT
WRITE(6,9901) IVEC
DO IALPH=1,NBAS,4
NREST=MIN(4,NBAS-IALPH+1)
DO III=1,NREST
IF (ABS(CI(IALPH-1+III,IVEC)).LT.1.D-9) THEN
CIOUT(III)=0.D0
ELSE
CIOUT(III)=CI(IALPH-1+III,IVEC)
END IF
END DO
WRITE(6,9902) (IALPH+III-1,CIOUT(III),III=1,NREST)
END DO
END DO
9901 FORMAT(' Function :',I5)
C 9902 FORMAT(5(I5,F10.5))
9902 FORMAT(4(I4,E14.6))
RETURN
END
C
SUBROUTINE RDMAT(NAME,ARRAY)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /UNITS/ IUNITR,IUNITX,IUNITO
C.. INCLUDE 'common_units.h'
DIMENSION ARRAY(NBASM,NBASM)
CHARACTER*7 NAME
C
DO I=1,NBAS
DO J=1,NBAS
ARRAY(I,J)=0.D0
END DO
END DO
C
OPEN(UNIT=IUNITR,FILE=NAME,STATUS='OLD',FORM='FORMATTED',ERR=8101)
100 CONTINUE
READ(IUNITR,*,IOSTAT=KK) I,J,XDUM
IF (KK.NE.0) GO TO 101
ARRAY(I,J)=XDUM
ARRAY(J,I)=XDUM
GO TO 100
101 CONTINUE
CLOSE(IUNITR)
RETURN
8101 CONTINUE
WRITE(6,*) ' ERROR WHILE OPENING FILE <',NAME,'>!! ERROR EXIT'
STOP ' ERROR IN RDMAT'
END
C
SUBROUTINE RDORS
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
C we have the "old" best density and the new density
C the density is PODA
COMMON /ODACOM/ PODA(NBASM,NBASM),XIODA(NBASM,NBASM),EMODA,EKODA
$ ,EHODA,EMPODA,LODAU
LOGICAL LODAU
C.. INCLUDE 'common_oda.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
CHARACTER*4 KEYW(2),STR3
CHARACTER*6 KEYOPT(30),STR6
CHARACTER*80 LINE
DATA KEYW /'*ORS','*END'/
DATA KEYOPT /'ECONVE','FCONVE','TEST ','CANONI','MAXITE',
- 'PIPEKM','LAMBDA','FREEZE','FMIXIN','REMOVE',
- 'WFN ','COREHA','MOLDEN','BOYS ','QMC ',
- 'COMPLE','EXTERN','MULTIP','ODA ','DIPOL ',
- 'ORIGIN','6 D ','REORDE','EXCHAN','ICSCF ',
- 'XXXXXX','XXXXXX','XXXXXX','XXXXXX','XXXXXX'/
C
C DEFAULTS
LODA=.FALSE.
LDIPOL=.FALSE.
LODAU=.TRUE.
LWFN=.FALSE.
LTST =.FALSE.
LECONV=.TRUE.
LCAN=.FALSE.
LPIPM=.FALSE.
LAMBDA=.FALSE.
NFRZ =0
XCMIX=0.30
LREMO=.FALSE.
LCOREH=.FALSE.
LMOLD=.FALSE.
LBOYS=.FALSE.
LQMC=.FALSE.
LCICOM=.FALSE.
LEXTER=.FALSE.
LMULTP=.FALSE.
LREORD=.FALSE.
L6D=.FALSE.
LICSCF=.FALSE.
XFAC=0.D0
DO I=1,3
ORIGIN(I)=0.D0
END DO
C
C structure of input:
C keyword for each program
IOINP=83
OPEN(UNIT=IOINP,FILE='INPUT.ORS',ERR=2217,FORM='FORMATTED',
- STATUS='OLD')
C first, look for keyword 'ORS'
1 CONTINUE
READ(IOINP,'(A4)',END=2,ERR=921) STR3
IF (STR3.EQ.KEYW(1)) THEN
C look for Keyoptions FCONVE, ECONVE, NITMAX, CANONI, TEST
11 CONTINUE
READ(IOINP,'(A6)',END=2,ERR=921) STR6
WRITE(6,*) ' read STR6 =|',STR6,'|'
IF (STR6(1:4).EQ.KEYW(2)) GO TO 2
IF (STR6.EQ.KEYOPT(1)) THEN
LECONV=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(2)) THEN
LECONV=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(3)) THEN
LTST=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(4)) THEN
LCAN=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(5)) THEN
READ(IOINP,*) NITMAX
ELSE IF (STR6.EQ.KEYOPT(6)) THEN
LPIPM=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(7)) THEN
C
C this lambda is an idea of Michel: we scale all bielectronic
c integrals by a common factor; lambda -> 0 gives the monoelectronic
c energy as total energy, and may be calculated by hand: it is the
c energy of the Bohr model -Z^2/n^2 in an infinite basis
c
LAMBDA=.TRUE.
READ(IOINP,*) XLAMB
ELSE IF (STR6.EQ.KEYOPT(8)) THEN
C
C freezing of orbitals in the Pipek-Mezey localization scheme
C
READ(IOINP,*) NFRZ
READ(IOINP,*) (IFRZ(I),I=1,NFRZ)
ELSE IF (STR6.EQ.KEYOPT(9)) THEN
READ(IOINP,*) XFMIX
IF (XFMIX.GT.1.D0) XFMIX=XFMIX*0.01D0
ELSE IF (STR6.EQ.KEYOPT(10)) THEN
LREMO=.TRUE.
READ(IOINP,*) NREMO
READ(IOINP,*) (IIREMO(I),I=1,NREMO)
ELSE IF (STR6.EQ.KEYOPT(11)) THEN
LWFN=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(12)) THEN
LCOREH=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(13)) THEN
LMOLD=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(14)) THEN
LBOYS=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(15)) THEN
LQMC=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(16)) THEN
LCICOM=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(17)) THEN
LEXTER=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(18)) THEN
LMULTP=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(19)) THEN
LODA=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(20)) THEN
LDIPOL=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(21)) THEN
READ(IOINP,*) (ORIGIN(I),I=1,3)
C L6D
ELSE IF (STR6.EQ.KEYOPT(22)) THEN
L6D=.TRUE.
LMULTP=.FALSE.
LPIPM=.FALSE.
LMOLD=.FALSE.
LQMC=.FALSE.
LBOYS=.FALSE.
LEXTER=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(23)) THEN
LREORD=.TRUE.
READ(IOINP,*) NATOM2
READ(IOINP,*) (IATSRT(I),I=1,NATOM2)
ELSE IF (STR6.EQ.KEYOPT(24)) THEN
C weight of exchange in density functional calculation
READ(IOINP,*) XFAC
ELSE IF (STR6.EQ.KEYOPT(25)) THEN
C ICSCF
LICSCF=.TRUE.
ELSE
WRITE(6,*)
WRITE(6,*) ' UNKNOWN OPTION: |',STR6,'|'
WRITE(6,*)
WRITE(6,*) ' POSSIBLE OPTIONS IN THE BLOCK *ORS ... *END ARE:'
WRITE(6,*)
WRITE(6,*) ' ',KEYOPT(1)
WRITE(6,*) ' ',KEYOPT(2)
WRITE(6,*) ' ',KEYOPT(3)
WRITE(6,*) ' ',KEYOPT(4)
WRITE(6,*) ' ',KEYOPT(5)
WRITE(6,*) ' ',KEYOPT(6)
WRITE(6,*) ' ',KEYOPT(7)
WRITE(6,*) ' ',KEYOPT(8)
WRITE(6,*) ' ',KEYOPT(9)
WRITE(6,*) ' ',KEYOPT(10)
WRITE(6,*) ' ',KEYOPT(11)
WRITE(6,*) ' ',KEYOPT(12)
WRITE(6,*) ' ',KEYOPT(13)
WRITE(6,*) ' ',KEYOPT(14)
WRITE(6,*) ' ',KEYOPT(15)
WRITE(6,*) ' ',KEYOPT(16)
WRITE(6,*) ' ',KEYOPT(17)
WRITE(6,*) ' ',KEYOPT(18)
WRITE(6,*) ' ',KEYOPT(19)
WRITE(6,*) ' ',KEYOPT(20)
WRITE(6,*) ' ',KEYOPT(21)
WRITE(6,*) ' ',KEYOPT(22)
WRITE(6,*) ' ',KEYOPT(23)
WRITE(6,*) ' ',KEYOPT(24)
WRITE(6,*) ' ',KEYOPT(25)
STOP ' CHOOSE CORRECT OPTION '
END IF
GO TO 11
END IF
GO TO 1
2 CONTINUE
CLOSE(IOINP)
RETURN
921 CONTINUE
CLOSE(IOINP)
STOP ' ERROR IN INPUT'
920 CONTINUE
CLOSE(IOINP)
STOP ' YOU SHOULD TERMINATE THE BLOCK ORT BY <*END> '
922 CONTINUE
CLOSE(IOINP)
STOP ' YOU SHOULD TERMINATE THE BLOCK SCI BY <*END> '
2217 CONTINUE
WRITE(6,*) ' NO FILE , USING THE DEFAULT VALUES'
RETURN
END
C
SUBROUTINE VSRTO(CI,IOCC,IOCCS)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
DIMENSION CI(NBASM,NBASM),IOCC(NBAS),IOCCS(NBAS)
C
IOCCZ=0
1 CONTINUE
IOCCZ=IOCCZ+1
IF (IOCCZ.EQ.NBAS+1) RETURN
IF (IOCC(IOCCZ).NE.IOCCS(IOCCZ)) THEN
IVIRTZ=IOCCZ
2 CONTINUE
IVIRTZ=IVIRTZ+1
IF (IVIRTZ.EQ.NBAS+1) RETURN
IF (IOCC(IVIRTZ).EQ.IOCCS(IVIRTZ)) THEN
IDUM=IOCC(IOCCZ)
IOCC(IOCCZ)=IOCC(IVIRTZ)
IOCC(IVIRTZ)=IDUM
IDUM=IOCCS(IOCCZ)
IOCCS(IOCCZ)=IOCCS(IVIRTZ)
IOCCS(IVIRTZ)=IDUM
DO J=1,NBAS
CDUM=CI(J,IOCCZ)
CI(J,IOCCZ)=CI(J,IVIRTZ)
CI(J,IVIRTZ)=CDUM
END DO
GO TO 1
END IF
GO TO 2
END IF
GO TO 1
END
SUBROUTINE WVEC(IORTYP)
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /UNITS/ IUNITR,IUNITX,IUNITO
C.. INCLUDE 'common_units.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
C
C output to a file
C
WRITE(6,*)
C
C localized orbitals by ORTHO
IF (IORTYP.EQ.1) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
WRITE(6,*) ' we assume ORTHO - localized orbitals '
IF (L6D) THEN
OPEN(UNIT=IUNITX,FILE='VECTOR.LOC.6D',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
ELSE
OPEN(UNIT=IUNITX,FILE='VECTOR.LOC',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
END IF
C canonical orbitals
ELSE IF (IORTYP.EQ.2) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
WRITE(6,*) ' we assume canonical orbitals '
IF (L6D) THEN
OPEN(UNIT=IUNITX,FILE='VECTOR.CAN.6D',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
ELSE
OPEN(UNIT=IUNITX,FILE='VECTOR.CAN',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
END IF
C Pipek-Mezey localized orbitals
ELSE IF (IORTYP.EQ.3) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
WRITE(6,*) ' we assume Pipek-Mezey localized orbitals '
OPEN(UNIT=IUNITX,FILE='VECTOR.PM',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
C Boys localized orbitals
ELSE IF (IORTYP.EQ.4) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
WRITE(6,*) ' we assume Boys localized orbitals '
OPEN(UNIT=IUNITX,FILE='VECTOR.BOYS',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
C temporary vector, for a crashed restart
ELSE IF (IORTYP.EQ.5) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
IF (L6D) THEN
OPEN(UNIT=IUNITX,FILE='VECTOR.TMP.6D',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
ELSE
OPEN(UNIT=IUNITX,FILE='VECTOR.TMP',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
END IF
C ICSCF orbitals
ELSE IF (IORTYP.EQ.6) THEN
WRITE(6,*) ' WRITING VECTOR TO FILE '
IF (L6D) THEN
OPEN(UNIT=IUNITX,FILE='VECTOR.ICSCF.6D',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
ELSE
OPEN(UNIT=IUNITX,FILE='VECTOR.ICSCF',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
END IF
ELSE
WRITE(6,*) ' demanded orbitals type not available ',IORTYP
WRITE(6,*) ' WRITING VECTOR TO FILE '
OPEN(UNIT=IUNITX,FILE='VECTOR.UNKNOWN',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
END IF
WRITE(6,*)
DO J=1,NBAS
WRITE(IUNITX,'(4E20.12)') (CI(K,J),K=1,NBAS)
WRITE(IUNITX,*)
END DO
WRITE(IUNITX,'(30I2)') (IOCC(I),I=1,NBAS)
WRITE(IUNITX,'(30I2)') (IOCCS(I),I=1,NBAS)
CLOSE(IUNITX)
C
RETURN
END
C
SUBROUTINE GENCAN(F,CI,ORBEN,NBAS)
INCLUDE 'param.h'
DIMENSION F(NBASM,NBASM),CI(NBASM,NBASM),ORBEN(NBASM)
DIMENSION CIN(NBASM,NBASM),CDUM1(NBASM),CDUM2(NBASM)
DIMENSION FDUM(NBASM,NBASM)
C
C the Fock matrix comes in in orthogonal molecular orbitals
C we give back the new orbitals and the new diagonal elements
C the Fock matrix is preserved
C
C here we diagonalize the Fock matrix, and store the new orbitals
C
DO I=1,NBAS
DO J=1,NBAS
FDUM(I,J)=F(I,J)
END DO
END DO
C
IONE=1
CALL RS(NBASM,NBAS,FDUM,ORBEN,IONE,CIN,CDUM1,CDUM2,IERR)
C
WRITE(6,*)
WRITE(6,*) ' ORBITAL ENERGIES OF THE CANONICAL ORBITALS'
WRITE(6,*)
WRITE(6,'(4(I4,E14.5))') (I,ORBEN(I),I=1,NBAS)
WRITE(6,*)
C
C well, it should come down to a simple matrix multiplication ...
C
DO IALPH=1,NBAS
DO I=1,NBAS
CDUM=0.D0
DO J=1,NBAS
C the J components of eigenvector I
CDUM=CDUM+CIN(J,I)*CI(IALPH,J)
END DO
FDUM(IALPH,I)=CDUM
END DO
END DO
C
DO IALPH=1,NBAS
DO I=1,NBAS
CI(IALPH,I)=FDUM(IALPH,I)
END DO
END DO
C
RETURN
END
C
SUBROUTINE PRECAN
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
C we have to obtain an orthonormal set of orbitals, and a start density
C matrix
IF (LCOREH) THEN
DO I=1,NBAS
DO J=I+1,NBAS
CI(I,J)=0.D0
CI(J,I)=0.D0
END DO
CI(I,I)=1.0D0
END DO
END IF
C
CALL SHALF(CI,NBAS,NBAS)
DO I=1,NBAS
DO J=1,NBAS
P(I,J)=0.0D0
END DO
END DO
C
IF (.NOT.LCOREH) THEN
C the density matrix
DO I=1,NBAS
DO IORB=1,NOCC
CII=2.D0*CI(I,IORB)
DO J=1,NBAS
P(I,J)=P(I,J)+CII*CI(J,IORB)
END DO
END DO
END DO
END IF
C
RETURN
END
C
SUBROUTINE PIPEKM
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
DIMENSION IKEEP(NBASM)
LOGICAL LSWEEP
DIMENSION HDUM(NBASM,NBASM)
C
C Pipek-Mezey localization, i.e. the occupied orbitals are localized,
C and as virtual orbitals the atomic orbitals are taken, and
C orthogonalized to the occupied and among themselves
C
WRITE(6,*)
WRITE(6,*) ' PIPEK - MEZEY LOCALIZATION : '
WRITE(6,*)
C
C orbital freezing: we put in first places all orbitals to be frozen
C
IF (NFRZ.NE.0) THEN
WRITE(6,*) ' we will freeze ',NFRZ,' orbitals'
WRITE(6,*)
IP1=0
IP2=1
1 CONTINUE
IF (IFRZ(IP2).EQ.1) THEN
IP1=IP1+1
IF (IP1.NE.IP2) THEN
DO I=1,NBAS
CDUM=CI(I,IP2)
CI(I,IP2)=CI(I,IP1)
CI(I,IP1)=CDUM
END DO
IDUM=IFRZ(IP2)
IFRZ(IP2)=IFRZ(IP1)
IFRZ(IP1)=IDUM
END IF
IP2=IP2+1
IF (IP2.LE.NBAS) GO TO 1
END IF
END IF
C
IF (NOCC-NFRZ.LE.1) THEN
WRITE(6,*) ' FOR NOCC=1 THERE IS NOTHING TO OPTIMIZE '
WRITE(6,*)
RETURN
END IF
CALL MEASUR(D,CI,NOCC)
WRITE(6,*) ' BEFORE: LOCALIZATION CRITERION D=',D
ITER=0
200 CONTINUE
ITER=ITER+1
CALL SWEEP(CI(1,NFRZ+1),NOCC-NFRZ,NBAS,NATOM,LSWEEP)
DOLD=D
CALL MEASUR(D,CI,NOCC)
c$$$ WRITE(6,*) ' LOCALIZATION CRITERION D=',D,' LSWEEP = ',LSWEEP
c$$$ IF (ABS(D-DOLD).LT.1.D-9) GO TO 201
IF (.NOT.LSWEEP) GO TO 201
GO TO 200
201 CONTINUE
WRITE(6,*) ' PIPEK-MEZEY PROCEDURE TOOK ',ITER
$ ,' SWEEPS TO CONVERGE'
WRITE(6,*) ' NOW: LOCALIZATION CRITERION D=',D
WRITE(6,*)
DO I=1,NBAS
DO J=1,NBAS
HDUM(I,J)=HAM(I,J)+XKIN(I,J)
END DO
END DO
CALL TRANS(HDUM,CI,NOCC,NBAS)
C reorder the vectors according to the orbital energy
DO I=1,NOCC
ORBEN(I)=HDUM(I,I)
END DO
C
C now we sort quite a lot at a time: ORBEN, CI, and F,
C by BUBBLESORT, despite the lecture of the Numerical Recipes
C
CALL BUBBLE(NOCC,NBAS,ORBEN,HDUM,CI)
C
WRITE(6,*) ' DIAGONAL ELEMENTS OF HCORE '
WRITE(6,'(4(I4,F12.6))') (I,HDUM(I,I),I=1,NOCC)
WRITE(6,*)
CALL CONVIR(CI)
C
RETURN
END
C
SUBROUTINE MEASUR(D,CI,NVECT)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
DIMENSION CI(NBASM,NVECT)
C
D=0.D0
XDUM=0.D0
DO I=1,NVECT
DI=0.D0
C construct the density matrix for orbital I
DO IAT=1,NATOM
QIA=0.D0
DO IALPH=ISTRT(IAT),ISTRT(IAT)+NABAS(IAT)
CIA=CI(IALPH,I)
DO IBETA=1,NBAS
PAB=CIA*CI(IBETA,I)
QIA=QIA+PAB*SAO(IBETA,IALPH)
END DO
END DO
DI=DI+QIA*QIA
XDUM=XDUM+QIA
END DO
D=D+DI
END DO
XDUM=2.D0*XDUM
C
C we should have XDUM = number of electrons
IF (ABS(2*NVECT-XDUM).GT.1.D-3) THEN
WRITE(6,*)
$ ' MEASUR: COUNTING ELECTRONS: SHOULD ',2*NVECT,' HAVE ',XDU
$M
C folded 1 (fixf $Revision: 1.3 $)
WRITE(6,*)
END IF
D=0.25D0*DBLE(2*NVECT)/D
RETURN
END
C
SUBROUTINE SWEEP(CI,NOCC,NBAS,NATOM,LSWEEP)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
DIMENSION SPS(NATMX),SPT(NATMX),TPT(NATMX),CI(NBASM,NBASM)
LOGICAL LSWEEP
C
SUMAL=0.D0
DO IS=1,NOCC
DO IT=IS+1,NOCC
C
C with every calculated angle orbital S will be changed as well
C we have to recalculate SPS again and again
C
DO IAT=1,NATOM
SSS=0.D0
SST=0.D0
STT=0.D0
DO IALPH=ISTRT(IAT),ISTRT(IAT)+NABAS(IAT)
CIAS=CI(IALPH,IS)
CIAT=CI(IALPH,IT)
DO IBETA=1,NBAS
CIBS=CI(IBETA,IS)
CIBT=CI(IBETA,IT)
SAOAB=SAO(IALPH,IBETA)
SSS=SSS + CIAS*CIBS *SAOAB
STT=STT + CIAT*CIBT *SAOAB
SST=SST + (CIAS*CIBT+CIAT*CIBS)*SAOAB
END DO
END DO
SPS(IAT)=SSS
SPT(IAT)=SST*0.5D0
TPT(IAT)=STT
END DO
AST=0.D0
BST=0.D0
DO IAT=1,NATOM
AST=AST+SPT(IAT)*SPT(IAT)-0.25D0*(SPS(IAT)-TPT(IAT))*(SPS(IAT)
$ -TPT(IAT))
BST=BST+SPT(IAT)*(SPS(IAT)-TPT(IAT))
END DO
S4A= BST/SQRT(AST*AST+BST*BST)
C4A=-AST/SQRT(AST*AST+BST*BST)
FOURA=ASIN(ABS(S4A))
PI=2.D0*ACOS(0.D0)
C which quadrant?
IF (S4A.GT.0) THEN
C between 0 and pi/2: do nothing
IF (C4A.LE.0) THEN
C between pi/2 and pi
FOURA=PI-FOURA
END IF
ELSE
IF (C4A.LT.0) THEN
C between pi and 3/2 pi
FOURA=FOURA+PI
ELSE
C between 3/2 pi and 2 pi
FOURA=2.D0*PI-FOURA
END IF
END IF
ALPHA=0.25D0*FOURA
ALPI=4.D0*ALPHA/PI
SUMAL=SUMAL+(NINT(ALPI)-ALPI)*(NINT(ALPI)-ALPI)
C
C there is this story about Gamma and Alpha in the paper of Pipek and
c Mezey; for BST=0 we obtain alpha=Pi/4, so it is gamma=alpha for
c localizing, thus minimizing D or maximizing P=D^-1
C
GAMMA=ALPHA
C
C now we have to mix the two orbitals
SINA=SIN(GAMMA)
COSA=COS(GAMMA)
DO IALPH=1,NBAS
CSA=CI(IALPH,IS)
CTA=CI(IALPH,IT)
CSSA=COSA*CSA+SINA*CTA
CTTA=COSA*CTA-SINA*CSA
CI(IALPH,IS)=CSSA
CI(IALPH,IT)=CTTA
END DO
C that's it
END DO
END DO
SUMAL=SQRT(SUMAL)/DBLE(NOCC*(NOCC-1)/2)
C if we need another sweep, we leave LSWEEP as true
LSWEEP=.TRUE.
C IF (SUMAL.LE.1.D-7) LSWEEP=.FALSE.
WRITE(6,*) ' SWEEP: SUMAL = ',SUMAL
IF (SUMAL.LE.1.D-9) LSWEEP=.FALSE.
C
RETURN
END
C
SUBROUTINE MOLDEN
C
C we try to create a MOLDEN file
C
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
DIMENSION IORD(NBASM)
CHARACTER*6 CHEXT
CHARACTER*2 SYMBAT(0:92)
DIMENSION NSHLT(NLMAX,NATMX),IZENTR(NSHLMX)
CHARACTER*1 CST(10),CBUCH(0:9)
CHARACTER*3 STRG(NBASM)
DATA CST /'S','P','D','F','G','H','I','J','K','L'/
DATA CBUCH /'0','1','2','3','4','5','6','7','8','9'/
DATA SYMBAT/ 'XX',' H','HE','LI','BE',' B',' C',' N',
1' O',' F','NE','NA','MG','AL','SI',' P',' S','CL','AR',
2' K','CA','SC','TI',' V','CR','MN','FE','CO','NI','CU',
3'ZN','GA','GE','AS','SE','BR','KR','RB','SR',' Y','ZR',
4'NB','MO','TC','RU','RH','PD','AG','CD','IN','SN','SB',
5'TE',' I','XE','CS','BA','LA','CE','PR','ND','PM','SM',
6'EU','GD','TB','DY','HO','ER','TM','YB','LU','HF','TA',
7' W','RE','OS','IR','PT','AU','HG','TL','PB','BI','PO',
8'AT','RN','FR','RA','AC','TH','PA',' U'/
WRITE(6,*)
WRITE(6,*) '-------------------------------------------------'
WRITE(6,*) '--- M O L D E N ----------- M O L D E N ---------'
WRITE(6,*) '-------------------------------------------------'
C
C
IOUQMC=57
OPEN(UNIT=IOUQMC,FILE='ORS.MOLDEN',FORM='FORMATTED',STATUS='UNKNOW
$N')
C folded 1 (fixf $Revision: 1.3 $)
CLOSE(IOUQMC,STATUS='DELETE')
OPEN(UNIT=IOUQMC,FILE='ORS.MOLDEN',FORM='FORMATTED',STATUS='NEW')
WRITE(IOUQMC,8001)
8001 FORMAT('[Molden Format]',/,'[Atoms] Angs')
C
IUNITR=76
C
C we read the file 2 times: first the atoms, then the basis set
C
OPEN(UNIT=IUNITR,FILE='SYSTEM.ORTHO',STATUS='OLD',
- FORM='FORMATTED')
C
WRITE(6,*)
C
C the first read, concerning atoms
C
READ(IUNITR,*) NATOM
DO IAT=1,NATOM
READ(IUNITR,*) NUMAT,NSHLAT,(POS(J,IAT),J=1,3)
WRITE(IOUQMC,8002) SYMBAT(NUMAT),IAT,NUMAT,(POS(J,IAT)*ANTOAU,
- J=1,3)
8002 FORMAT(A2,2I5,F21.10,2F20.10)
C
DO ISH=1,NSHLAT
READ(IUNITR,*) IDUM,NNPRIM
DO III=1,NNPRIM
READ(IUNITR,*) XDUM
END DO
END DO
END DO
WRITE(IOUQMC,8003)
8003 FORMAT('[5D]',/,'[GTO]')
C
C the second read of SYSTEM.ORTHO
C
REWIND(IUNITR)
NSHL=0
NBASY=1
LMAX=0
IPRIM=0
NPRIMT=0
IBAS=0
READ(IUNITR,*) NATOM
DO IAT=1,NATOM
DO ITYPE=1,NLMAX
NSHLT(ITYPE,IAT)=0
END DO
READ(IUNITR,*) IDUM,NSHLAT
WRITE(IOUQMC,8013) IAT
8013 FORMAT(I3,' 0')
C
DO ISH=1,NSHLAT
NSHL=NSHL+1
IZENTR(NSHL)=IAT
READ(IUNITR,*) ITYPE,NPRIM(NSHL)
WRITE(IOUQMC,8004) CST(ITYPE+1),NPRIM(NSHL)
8004 FORMAT(1X,A1,I5,' 1.00')
NPRIMT=NPRIMT+(ITYPE+ITYPE+1)*NPRIM(NSHL)
IL(NSHL)=ITYPE
LMAX=MAX(ITYPE+1,LMAX)
NSHLT(ITYPE+1,IAT)=NSHLT(ITYPE+1,IAT)+1
NBASY=NBASY+ITYPE+ITYPE+1
IPST=IPRIM+1
DO III=1,NPRIM(NSHL)
IPRIM=IPRIM+1
READ(IUNITR,*) EXX(IPRIM),COEFF(IPRIM)
END DO
C
C we have to reorder the indices of the primitives for having the MOLDEN
C ordering of m values
C
C we, i.e. DALTON:
C s px py pz d-2 d-1 d0 d1 d2 f-3 f-2 f-1 f0 f1 f2 f3
C
C MOLDEN likes the order
C s px py pz d0 d1 d-1 d2 d-2 f0 f1 f-1 f2 f-2 f3 f-3
C
C 1 | 1 2 3 | 1 2 3 4 5 | 1 2 3 4 5 6 7
C 1 | 3 1 2 | 3 4 2 5 1 | 4 5 3 6 2 7 1
C
IF (ITYPE.EQ.0) THEN
C s functions
IBAS=IBAS+1
IORD(IBAS)=IBAS
STRG(IBAS)='S '
ELSE IF (ITYPE.EQ.1) THEN
C p functions
IORD(IBAS+1)=IBAS+1
STRG(IBAS+1)='PX '
IORD(IBAS+2)=IBAS+2
STRG(IBAS+2)='PY '
IORD(IBAS+3)=IBAS+3
STRG(IBAS+3)='PZ '
IBAS=IBAS+3
ELSE IF (ITYPE.EQ.2) THEN
C d functions
IORD(IBAS+5)=IBAS+1
STRG(IBAS+5)='D-2'
IORD(IBAS+3)=IBAS+2
STRG(IBAS+3)='D-1'
IORD(IBAS+1)=IBAS+3
STRG(IBAS+1)='D 0'
IORD(IBAS+2)=IBAS+4
STRG(IBAS+2)='D 1'
IORD(IBAS+4)=IBAS+5
STRG(IBAS+4)='D 2'
IBAS=IBAS+5
ELSE IF (ITYPE.EQ.3) THEN
C f functions
IORD(IBAS+7)=IBAS+1
STRG(IBAS+7)='F-3'
IORD(IBAS+5)=IBAS+2
STRG(IBAS+5)='F-2'
IORD(IBAS+3)=IBAS+3
STRG(IBAS+3)='F-1'
IORD(IBAS+1)=IBAS+4
STRG(IBAS+1)='F 0'
IORD(IBAS+2)=IBAS+5
STRG(IBAS+2)='F 1'
IORD(IBAS+4)=IBAS+6
STRG(IBAS+4)='F 2'
IORD(IBAS+6)=IBAS+7
STRG(IBAS+6)='F 3'
IBAS=IBAS+7
ELSE
STOP ' no functions beyond l=3 installed yet '
END IF
C
CALL NORMBF(NSHL,IPST)
DO III=0,NPRIM(NSHL)-1
WRITE(IOUQMC,'(2F18.10)') EXX(IPST+III),COEFF(IPST+III)
END DO
END DO
WRITE(IOUQMC,*)
END DO
C
CLOSE(IUNITR)
WRITE(IOUQMC,*)
WRITE(IOUQMC,8005)
8005 FORMAT('[MO]')
C
WRITE(6,*) ' NBAS =',NBAS
WRITE(6,*) ' NOCC =',NOCC
WRITE(6,*) ' LMAX =',LMAX
WRITE(6,*) ' NSHL =',NSHL
WRITE(6,*) ' NPRIM =',(NPRIM(I),I=1,NSHL)
WRITE(6,*) ' NPRIMT =',NPRIMT
C
C write the MOs
C
DO I=1,NBAS
WRITE(IOUQMC,8006) ORBEN(I),IOCC(I)
8006 FORMAT(' Ene=',F13.4,/,' Spin= Alpha',/,' Occup=',I4,'.000000')
C well,z not correct neither
c$$$ WRITE(IOUQMC,'(I4,9H ,A3,F15.6)') (IALPH,STRG(IALPH)
c$$$ $ ,CI(IORD(IALPH),I),IALPH=1,NBAS)
WRITE(IOUQMC,'(I4,F15.6,3H ,A3)') (IALPH,CI(IORD(IALPH),I)
$ ,STRG(IALPH),IALPH=1,NBAS)
END DO
C that's it
CLOSE(IOUQMC)
RETURN
END
C
C
SUBROUTINE SWEEPB(HERMAT,CI,NOCC,NBAS,LSWEEP)
INCLUDE 'param.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C the common things like NFRZ or NREMO we take from the Pipek-Mezey
C localization
C.. INCLUDE 'common_boys.h'
DIMENSION CI(NBASM,NBASM)
DIMENSION HERMAT(NBASM,NBASM)
LOGICAL LSWEEP
C
A0=0.D0
A1=1.D0
A2=2.D0
A10=10.D0
AQU=0.25D0
THSW=A0
DO 110 I=2,NOCC
C
IM1=I-1
IP1=I+1
DO 100 J=1,IM1
C
C LOOP OVER I,J ROTATIONS
C MAXIMIZE SUM(K) **2 + **2
C LIKE IN JACOBI
C
C GET AX , B DETERMINING THE TAN OF THE ROTATION ANGLE
C
AX=A0
B=A0
DO K=1,3
AUX1=DIPOL(J,J,K)-DIPOL(I,I,K)
AUX2=DIPOL(I,J,K)
AX=AX+AUX1*AUX2
B=B+AQU*AUX1**2-AUX2**2
END DO
C
C ROTATION ANGLE X0
C
IF (DABS(AX)+DABS(B).LT.1.D-12) GO TO 100
X0=AQU*DATAN2(AX,B)
c$$$
c$$$ WRITE(6,*) ' I, J, AX, B, X0',I,J,AX,B,X0
c$$$ WRITE(6,*) DIPOL(I,I,1),DIPOL(J,J,1),DIPOL(I,J,1)
THSW=DMAX1(THSW,DABS(X0))
IF (DABS(X0).LT.1.D-13) GO TO 100
C
C ROTATION
C
IT=IT+1
C C=DCOS(X0)
S=DSIN(X0)
S2=S*S
C2=A1-S2
C=DSQRT(C2)
C
C update of the rotation matrix
C
DO 40 L=1,NOCC
AUX1=HERMAT(L,J)
HERMAT(L,J)=C*HERMAT(L,J)+S*HERMAT(L,I)
HERMAT(L,I)=C*HERMAT(L,I)-S*AUX1
40 CONTINUE
C
C update of the vector
C
DO IALPH=1,NBAS
CSA=CI(IALPH,I)
CTA=CI(IALPH,J)
CI(IALPH,I)=C*CSA-S*CTA
CI(IALPH,J)=C*CTA+S*CSA
END DO
C
JM1=J-1
C C2=C*C
C S2=S*S
CS=C*S
C2MS2=C2-S2
JP1=J+1
DO 95 K=1,3
AUX1=DIPOL(I,I,K)
AUX2=DIPOL(J,J,K)
DIPOL(I,I,K)=AUX1*C2+AUX2*S2-A2*DIPOL(I,J,K)*CS
DIPOL(J,J,K)=AUX1*S2+AUX2*C2+A2*DIPOL(I,J,K)*CS
DIPOL(I,J,K)=CS*(AUX1-AUX2)+DIPOL(I,J,K)*C2MS2
DIPOL(J,I,K)=DIPOL(I,J,K)
DO L=1,JM1
AUX1=C*DIPOL(I,L,K)-S*DIPOL(J,L,K)
DIPOL(J,L,K)=S*DIPOL(I,L,K)+C*DIPOL(J,L,K)
DIPOL(I,L,K)=AUX1
DIPOL(L,I,K)=DIPOL(I,L,K)
DIPOL(L,J,K)=DIPOL(J,L,K)
END DO
DO L=JP1,IM1
AUX1=C*DIPOL(I,L,K)-S*DIPOL(J,L,K)
DIPOL(J,L,K)=C*DIPOL(J,L,K)+S*DIPOL(I,L,K)
DIPOL(I,L,K)=AUX1
DIPOL(L,I,K)=DIPOL(I,L,K)
DIPOL(L,J,K)=DIPOL(J,L,K)
END DO
DO L=IP1,NOCC
AUX1=C*DIPOL(J,L,K)+S*DIPOL(I,L,K)
DIPOL(I,L,K)=C*DIPOL(I,L,K)-S*DIPOL(J,L,K)
DIPOL(J,L,K)=AUX1
DIPOL(L,I,K)=DIPOL(I,L,K)
DIPOL(L,J,K)=DIPOL(J,L,K)
END DO
95 CONTINUE
100 CONTINUE
110 CONTINUE
C
C CHECK CONVERGENCE
C
C if we need another sweep, we leave LSWEEP as true
LSWEEP=.TRUE.
WRITE(6,*) ' SWEEP: GAMMAX = ',THSW
IF (THSW.LE.1.D-9) LSWEEP=.FALSE.
C
RETURN
END
C
SUBROUTINE BOYS
C
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C the common things like NFRZ or NREMO we take from the Pipek-Mezey
C localization
C.. INCLUDE 'common_boys.h'
DIMENSION IKEEP(NBASM),DIPDIA(NBASM,3)
LOGICAL LSWEEP
DIMENSION HDUM(NBASM,NBASM),ORD(NBASM)
DIMENSION HERMAT(NBASM,NBASM)
C
C Boys localization, i.e. the occupied orbitals are localized,
C and as virtual orbitals the atomic orbitals are taken, and
C orthogonalized to the occupied and among themselves
C
C
C we follow the recipe of Ahlrichs from TURBOMOLE:
C condition for minimization of \sum_i**2
C is *(- = 0 for i.ne.j
C so we set up the matrix of conditions
C
WRITE(6,*)
WRITE(6,*) ' BOYS LOCALIZATION : '
WRITE(6,*)
NVIRT=NBAS-NOCC
C
C orbital freezing: we put in first places all orbitals to be frozen
C
IF (NFRZ.NE.0) THEN
WRITE(6,*) ' we will freeze ',NFRZ,' orbitals'
WRITE(6,*)
IP1=0
IP2=1
1 CONTINUE
IF (IFRZ(IP2).EQ.1) THEN
IP1=IP1+1
IF (IP1.NE.IP2) THEN
DO I=1,NBAS
CDUM=CI(I,IP2)
CI(I,IP2)=CI(I,IP1)
CI(I,IP1)=CDUM
END DO
IDUM=IFRZ(IP2)
IFRZ(IP2)=IFRZ(IP1)
IFRZ(IP1)=IDUM
END IF
IP2=IP2+1
IF (IP2.LE.NBAS) GO TO 1
END IF
END IF
C
IF (NOCC-NFRZ.LE.1) THEN
WRITE(6,*) ' FOR NOCC=1 THERE IS NOTHING TO OPTIMIZE '
WRITE(6,*)
RETURN
END IF
IPRI=1
CALL MEASUR(D,CI,NOCC)
C
C get the dipole matrix elements
C
CALL RDMAT('DIPOL_X',DIPOL(1,1,1))
CALL RDMAT('DIPOL_Y',DIPOL(1,1,2))
CALL RDMAT('DIPOL_Z',DIPOL(1,1,3))
CALL TRANS(DIPOL(1,1,1),CI(1,NFRZ+1),NOCC-NFRZ,NBAS)
CALL TRANS(DIPOL(1,1,2),CI(1,NFRZ+1),NOCC-NFRZ,NBAS)
CALL TRANS(DIPOL(1,1,3),CI(1,NFRZ+1),NOCC-NFRZ,NBAS)
DO I=1,NOCC-NFRZ
DO K=1,3
DIPDIA(I,K)=DIPOL(I,I,K)
END DO
END DO
C
C we evaluate the Boys measure:
CALL BOYSME(NOCC-NFRZ,DIPDIA)
C
ITER=0
200 CONTINUE
ITER=ITER+1
C
C eigentlich sollte eine Boys Lokalisierung nur eine andere Art eines
C SWEEP sein ....
C
CALL SWEEPB(HERMAT,CI(1,NFRZ+1),NOCC-NFRZ,NBAS,LSWEEP)
C
C CALL TSTUNI(HERMAT,NBASM,NOCC-NFRZ)
C
IF (.NOT.LSWEEP) GO TO 201
GO TO 200
201 CONTINUE
WRITE(6,*) ' BOYS PROCEDURE TOOK ',ITER,' SWEEPS TO CONVERGE'
WRITE(6,*)
C we evaluate the Boys measure:
DO I=1,NOCC-NFRZ
DO K=1,3
DIPDIA(I,K)=DIPOL(I,I,K)
END DO
END DO
CALL BOYSME(NOCC-NFRZ,DIPDIA)
C
D=0.D0
CALL MEASUR(D,CI,NOCC)
C
IF (NFRZ.NE.0) THEN
C
C we have to take into account the frozen orbitals
C
DO I=NOCC,NFRZ,-1
DO J=NOCC,NFRZ,-1
HERMAT(I,J)=HERMAT(I-NFRZ,J-NFRZ)
HERMAT(I-NFRZ,J-NFRZ)=0.D0
END DO
END DO
DO I=1,NFRZ
HERMAT(I,I)=1.D0
END DO
END IF
DO I=1,NBAS
DO J=1,NBAS
HDUM(I,J)=HAM(I,J)+XKIN(I,J)
END DO
END DO
CALL TRANS(HDUM,CI,NOCC,NBAS)
C reorder the vectors according to the orbital energy
DO I=1,NOCC
ORD(I)=HDUM(I,I)
END DO
C
C now we sort quite a lot at a time: ORBEN, CI, and F,
C by BUBBLESORT, despite the lecture of the Numerical Recipes
C
CALL BUBBLE(NOCC,NBAS,ORD,HDUM,CI)
C
WRITE(6,*)
C
WRITE(6,*) ' DIAGONAL ELEMENTS OF HCORE '
WRITE(6,'(4(I4,F12.6))') (I,HDUM(I,I),I=1,NOCC)
WRITE(6,*)
C
C construct the virtual space
C
CALL CONVIR(CI)
C
RETURN
END
C
SUBROUTINE BOYSME(NORB,DIPDIA)
INCLUDE 'param.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C the common things like NFRZ or NREMO we take from the Pipek-Mezey
C localization
C.. INCLUDE 'common_boys.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
DIMENSION DIPDIA(NBASM,3)
C get the barycenter
XC1=0.D0
YC1=0.D0
ZC1=0.D0
DO I=1,NATOM
XC1=XC1+POS(1,I)
YC1=YC1+POS(2,I)
ZC1=ZC1+POS(3,I)
END DO
XC1=XC1/DBLE(NATOM)
YC1=YC1/DBLE(NATOM)
ZC1=ZC1/DBLE(NATOM)
C
WRITE(6,9771) NATOM,XC1,YC1,ZC1
9771 FORMAT(' barycenter of the molecule (',I2,' atoms) : ',3F10.4)
C
C print the list of orbital centers
C
WRITE(6,9909)
9909 FORMAT(/,' barycenters of the orbitals '
- ,//,5x,'No',10x,'x',15x,'y',15x,'z',/)
DO I=1,NORB
WRITE(6,9902) I,(DIPDIA(I,K),K=1,3)
9902 FORMAT(I6,3F16.8)
END DO
BBB=0.D0
DO I=1,NORB
DO J=I+1,NORB
DO K=1,3
BBB=BBB
- +(DIPDIA(I,K)-DIPDIA(J,K))*(DIPDIA(I,K)-DIPDIA(J,K))
END DO
END DO
END DO
C
WRITE(6,*)
WRITE(6,*) ' LOCALIZATION CRITERION =',BBB
WRITE(6,*)
C
RETURN
END
SUBROUTINE CONVIR(CI)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /MEZEY/ IFRZ(NBASM),ISTRT(NATMX),NABAS(NATMX)
$ ,NFRZ,IIREMO(NBASM),NREMO,ITMP(NBASM)
C.. INCLUDE 'common_mezey.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
DIMENSION CI(NBASM,NBASM)
LOGICAL LSWEEP
DIMENSION CIVIRT(NBASM,NBASM),CIDIFF(NBASM),CIDUM(NBASM,NBASM)
DIMENSION SDIAG(NBASM),CDUM1(NBASM),CDUM2(NBASM),SMO(NBASM,NBASM)
DIMENSION SPREAD(NBASM)
C
C DELCP
WRITE(6,*)
WRITE(6,*) ' CONSTRUCTING THE VIRTUAL SPACE '
WRITE(6,*)
C
C first of all we construct the NBAS basis functions without the
C occupied space
C
DO J=1,NBAS
IBETA=J
DO IALPH=1,NBAS
CIDIFF(IALPH)=0.D0
END DO
DO I=1,NOCC
C calculate the overlap
C CIVIRT is the unit matrix
SIJ=0.D0
DO IALPH=1,NBAS
SIJ=SIJ+CI(IALPH,I)*SAO(IALPH,IBETA)
END DO
C get the difference vector
DO IALPH=1,NBAS
CIDIFF(IALPH)=CIDIFF(IALPH)-SIJ*CI(IALPH,I)
END DO
END DO
C project
DO IALPH=1,NBAS
CIVIRT(IALPH,J)=CIDIFF(IALPH)
END DO
C add the diagonal, i.e. the orginal vector, IBETA=J
CIVIRT(IBETA,J)=CIVIRT(IBETA,J)+1.D0
END DO
C
C--SKIPCP
C this is the systematic, however delocalizing construction
C
C transform the overlap matrix
DO I=1,NBAS
DO J=1,NBAS
SMO(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(SMO,CIVIRT,NBAS,NBAS)
C
C diagonalize the overlap matrix
C
IONE=1
CALL RS(NBASM,NBAS,SMO,SDIAG,IONE,CIDUM,CDUM1,CDUM2,IERR)
C
WRITE(6,*)
WRITE(6,*) ' EIGENVALUES OF THE PROJECTED OVERLAP MATRIX'
WRITE(6,*)
WRITE(6,'(4(I4,F14.8))') (I,SDIAG(I),I=1,NBAS)
WRITE(6,*)
C
C the first NVIRT orbitals are the new virtual orbitals,
C there should be zero for the first NOCC eigenvalues
C
DO I=NOCC+1,NBAS
DO IALPH=1,NBAS
CIDIFF(IALPH)=0.D0
END DO
DO J=1,NBAS
CIJ=CIDUM(J,I)
DO IALPH=1,NBAS
CIDIFF(IALPH)=CIDIFF(IALPH)+CIVIRT(IALPH,J)*CIJ
END DO
END DO
DO IALPH=1,NBAS
CI(IALPH,I)=CIDIFF(IALPH)
END DO
END DO
C
C we need to normalize the new virtual orbitals for some reasons
CALL SNORM(CI(1,NOCC+1),NBAS-NOCC,NBAS)
C
CALL MEASUR(D,CI(1,NOCC+1),NBAS-NOCC)
WRITE(6,*) ' BEFORE: LOCALIZATION CRITERION D=',D
ITER=0
200 CONTINUE
ITER=ITER+1
CALL SWEEP(CI(1,NOCC+1),NBAS-NOCC,NBAS,NATOM,LSWEEP)
DOLD=D
CALL MEASUR(D,CI(1,NOCC+1),NBAS-NOCC)
c$$$ WRITE(6,*) ' LOCALIZATION CRITERION D=',D,' LSWEEP = ',LSWEEP
C
c$$$ IF (ABS(D-DOLD).LT.1.D-9) GO TO 201
IF (.NOT.LSWEEP) GO TO 201
GO TO 200
201 CONTINUE
WRITE(6,*) ' PIPEK-MEZEY PROCEDURE TOOK ',ITER
$ ,' SWEEPS TO CONVERGE'
WRITE(6,*) ' NOW: LOCALIZATION CRITERION D=',D
WRITE(6,*)
C
C//ENDCP
C
CP DO I=1,NBAS
CP DI=0.D0
C construct the density matrix for orbital I
CP DO IAT=1,NATOM
CP QIA=0.D0
CP DO IALPH=ISTRT(IAT),ISTRT(IAT)+NABAS(IAT)
CP CIA=CIVIRT(IALPH,I)
CP DO IBETA=1,NBAS
CP PAB=CIA*CIVIRT(IBETA,I)
CP QIA=QIA+PAB*SAO(IBETA,IALPH)
CP END DO
CP END DO
CP DI=DI+QIA*QIA
CP XDUM=XDUM+QIA
CP END DO
CP SPREAD(I)=DI
CP END DO
C
CP WRITE(6,*)
CP WRITE(6,*) ' SPREAD OF THE CONSTRUCTED ORBITALS '
CP WRITE(6,*)
CP WRITE(6,'(4(I4,F14.8))') (I,SPREAD(I),I=1,NBAS)
CP WRITE(6,*)
C
CP CALL BUBBLE(NBAS,NBAS,SPREAD,SMO,CIVIRT)
CP WRITE(6,*)
CP WRITE(6,*) ' SPREAD OF THE CONSTRUCTED ORBITALS '
CP WRITE(6,*)
CP WRITE(6,'(4(I4,F14.8))') (I,SPREAD(I),I=1,NBAS)
CP WRITE(6,*)
C
CP STOP
c$$$
c$$$C
c$$$C we may delete just the AO with the largest coefficient, one for each
c$$$C occupied MO
c$$$C
c$$$ DO I=1,NBAS
c$$$ IKEEP(I)=1
c$$$ END DO
c$$$
c$$$ WRITE(6,*) ' LREMO is set to ',LREMO
c$$$C
c$$$ IF (LREMO) THEN
c$$$C we remove predefined orbitals
c$$$ WRITE(6,*) ' removing orbitals by predefined list: '
c$$$ DO I=1,NOCC
c$$$ WRITE(6,9001) I,IIREMO(I)
c$$$ IKEEP(IIREMO(I))=0
c$$$ END DO
c$$$ ELSE
c$$$ DO IVECT=1,NOCC
c$$$ XL=0.D0
c$$$ DO IALPH=1,NBAS
c$$$ IF (ABS(CI(IALPH,IVECT)).GT.XL.AND.IKEEP(IALPH).EQ.1) THEN
c$$$ XL=ABS(CI(IALPH,IVECT))
c$$$ IREMO=IALPH
c$$$ END IF
c$$$ END DO
c$$$ WRITE(6,9001) IVECT,IREMO
c$$$ 9001 FORMAT(' OCC No',I4,': REMOVING AO No ',I4)
c$$$ IKEEP(IREMO)=0
c$$$ END DO
c$$$ END IF
c$$$
c$$$ DO IVECT=NOCC+1,NBAS
c$$$ DO IALPH=1,NBAS
c$$$ CI(IALPH,IVECT)=0.D0
c$$$ END DO
c$$$
c$$$ IO=1
c$$$ 100 CONTINUE
c$$$ IF (IKEEP(IO).EQ.1) THEN
c$$$ IKEEP(IO)=0
c$$$ GO TO 101
c$$$ END IF
c$$$ IO=IO+1
c$$$ GO TO 100
c$$$ 101 CONTINUE
c$$$ CI(IO,IVECT)=1.D0
c$$$C
c$$$ END DO
c$$$C
c$$$C we orthogonalize
c$$$c$$$ CALL SHALF(CI,NOCC,NBAS)
c$$$ WRITE(6,*)
c$$$ WRITE(6,*) ' ORTHOGONALIZING THE O-V block '
c$$$ CALL SFULL(CI,NOCC,NBAS)
c$$$ WRITE(6,*)
c$$$ WRITE(6,*) ' ORTHOGONALIZING THE V-V block '
c$$$ CALL SHALF(CI(1,NOCC+1),NBAS-NOCC,NBAS)
RETURN
END
C
SUBROUTINE OUTQMC
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
C
COMMON /POLYTA/ FAC(NLMAX+NLMAX),FACC(NLMAX)
C
LOGICAL LI,LJ,LK,LL
DIMENSION IVALPH(NBASM),IVBETA(NBASM)
C
DIMENSION OWBUF(NPRIMX*(NLMAX+NLMAX+1)),
- IWBUF(NPRIMX*(NLMAX+NLMAX+1)),XNBUF(NPRIMX*(NLMAX+NLMAX+1))
DIMENSION XNORM(NLMAX,-NLMAX:NLMAX),DNORM(NPRIMX*(NLMAX+NLMAX+1))
DIMENSION XCOEF(NLMAX,-NLMAX:NLMAX)
CHARACTER*16 CSTRI(NLMAX,-NLMAX:NLMAX)
CHARACTER*2 SYMBAT(0:92)
INTEGER*2 SYMLEN(0:92)
DIMENSION NSHLT(NLMAX,NATMX),COPRIM(NPRIMX),IZENTR(NSHLMX)
CHARACTER*1 CST(10),CBUCH(0:9)
DATA CST /'S','P','D','F','G','H','I','J','K','L'/
DATA CBUCH /'0','1','2','3','4','5','6','7','8','9'/
DATA SYMBAT/ 'XX','H ','HE','LI','BE','B ','C ','N ',
1'O ','F ','NE','NA','MG','AL','SI','P ','S ','CL','AR',
2'K ','CA','SC','TI','V ','CR','MN','FE','CO','NI','CU',
3'ZN','GA','GE','AS','SE','BR','KR','RB','SR','Y ','ZR',
4'NB','MO','TC','RU','RH','PD','AG','CD','IN','SN','SB',
5'TE','I ','XE','CS','BA','LA','CE','PR','ND','PM','SM',
6'EU','GD','TB','DY','HO','ER','TM','YB','LU','HF','TA',
7'W ','RE','OS','IR','PT','AU','HG','TL','PB','BI','PO',
8'AT','RN','FR','RA','AC','TH','PA','U '/
DATA SYMLEN/ 2,1,2,2,2,1,1,1,
1 1,1,2,2,2,2,2,1,1,2,2,
2 1,2,2,2,1,2,2,2,2,2,2,
3 2,2,2,2,2,2,2,2,2,1,2,
4 2,2,2,2,2,2,2,2,2,2,2,
5 2,1,2,2,2,2,2,2,2,2,2,
6 2,2,2,2,2,2,2,2,2,2,2,
7 1,2,2,2,2,2,2,2,2,2,2,
8 2,2,2,2,2,2,2,1/
CSTRI(1,0) = 'S 1 '
CSTRI(2,-1) = 'P Y '
CSTRI(2,1) = 'P X '
CSTRI(2,0) = 'P Z '
CSTRI(3, 2) = 'D XX-YY '
CSTRI(3, 1) = 'D XZ '
CSTRI(3, 0) = 'D 3ZZ-RR '
CSTRI(3,-1) = 'D YZ '
CSTRI(3,-2) = 'D XY '
CSTRI(4, 3) = 'F X(XX-3YY) '
CSTRI(4, 2) = 'F Z(XX-YY) '
CSTRI(4, 1) = 'F X(5ZZ-RR) '
CSTRI(4, 0) = 'F Z(5ZZ-RR) '
CSTRI(4,-1) = 'F Y(5ZZ-RR) '
CSTRI(4,-2) = 'F XYZ '
CSTRI(4,-3) = 'F Y(3XX-YY) '
CSTRI(5, 4)='G X(XX-3YY) '
CSTRI(5, 3)='G X(XX-3YY) '
CSTRI(5, 2)='G Z(XX-YY) '
CSTRI(5, 1)='G X(5ZZ-RR) '
CSTRI(5, 0)='G Z(5ZZ-RR) '
CSTRI(5,-1)='G Y(5ZZ-RR) '
CSTRI(5,-2)='G XYZ '
CSTRI(5,-3)='G Y(3XX-YY) '
CSTRI(5,-4)='G X(XX-3YY) '
XCOEF(1, 0)= 1.0D0
XCOEF(2,-1)= 1.0D0
XCOEF(2, 0)= 1.0D0
XCOEF(2, 1)= 1.0D0
XCOEF(3, 2)= 3.0D0
XCOEF(3, 1)= 3.0D0
XCOEF(3, 0)= 0.5D0
XCOEF(3,-1)= 3.0D0
XCOEF(3,-2)= 6.0D0
XCOEF(4, 3)=15.0D0
XCOEF(4, 2)=15.0D0
XCOEF(4, 1)= 1.5D0
XCOEF(4, 0)= 0.5D0
XCOEF(4,-1)= 1.5D0
XCOEF(4,-2)=30.0D0
XCOEF(4,-3)=15.0D0
WRITE(6,*)
WRITE(6,*) '-------------------------------------------------'
WRITE(6,*) '--- Q M C ----------- Q M C----------------------'
WRITE(6,*) '-------------------------------------------------'
C
C this piece of code is from ICMP.
C
C
C two files for the same information are opened:
C QMCDETLST_???? contains all relevant data and
C
C qmcin_???? contains the input file for the program of Michel
C
C the files are of the following structure:
C
C QMCDETLST_????: number of atoms, positions, basis, determiniants
C
C qmcin_????: number of atoms, positions, #determinants, Jastrows,
C basis, determinants
C
C so we determine already here the number of determinants in the
C expansion
C
C
IOUQMC=57
IINQMC=59
OPEN(UNIT=IOUQMC,FILE='QMCDETLST',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
CLOSE(IOUQMC,STATUS='DELETE')
OPEN(UNIT=IOUQMC,FILE='QMCDETLST',FORM='FORMATTED',STATUS
$ ='NEW')
OPEN(UNIT=IINQMC,FILE='qmcin',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
CLOSE(IINQMC,STATUS='DELETE')
OPEN(UNIT=IINQMC,FILE='qmcin',FORM='FORMATTED',STATUS
$ ='NEW')
WRITE(IOUQMC,*) ' Information written for a Hartree-Fock run'
WRITE(IOUQMC,*) ' ============================================'
WRITE(IOUQMC,*)
C
C we give the whole information to a file, i.e. basis set, position and
C kind of atoms ... for this purpose we have to reread the system input,
C normalize basis functions, break Gaussians down ....
C
CALL FFCAL
C
C all the basis set data is present from the beginning
C
C IOUQMC is a more general file with all necessary information
WRITE(IOUQMC,8001) NATOM
C IINQMC is the input file for the program of Michel Caffarel
WRITE(IINQMC,8101) NATOM
8001 FORMAT(' Number of atoms in the Molecule',/,I5)
8101 FORMAT('************************************',/
$ ,'Begin: Describe molecule',/
$ ,'************************************',/,'Number of nuclei',
$ /,I4,/,'Charge and position of nuclei')
WRITE(6,*)
C
NUMDET=1
NSHL=0
NBASY=0
LMAX=0
IPRIM=0
NPRIMT=0
DO IAT=1,NATOM
DO ITYPE=1,NLMAX
NSHLT(ITYPE,IAT)=0
END DO
NUMAT=NZ(IAT)
WRITE(IOUQMC,8002) SYMBAT(NUMAT),NUMAT,(POS(J,IAT),J=1,3)
WRITE(IINQMC,8102) NUMAT,(POS(J,IAT),J=1,3)
8002 FORMAT(' Atomic Symbol :',A,' nuclear charge '/,I5,/
$ ,' Position of this atom in a.u.',/,3(F20.12),/)
8102 FORMAT(I4,3F20.12)
C
WRITE(IOUQMC,8003)
8003 FORMAT(' Original Contracted Basis of that atom,'
$ ,' functions normalized (w.r.t. primtives)')
DO ISH=1,NSH(IAT)
NSHL=NSHL+1
IZENTR(NSHL)=IAT
ITYPE=IL(NSHL)
NPRIMT=NPRIMT+(ITYPE+ITYPE+1)*NPRIM(NSHL)
WRITE(IOUQMC,8004) ITYPE,NPRIM(NSHL),NSHL,CST(ITYPE+1)
8004 FORMAT(2I5,' ( contraction No',I4,' type ',A,' )')
LMAX=MAX(ITYPE+1,LMAX)
NSHLT(ITYPE+1,IAT)=NSHLT(ITYPE+1,IAT)+1
NBASY=NBASY+ITYPE+ITYPE+1
IPST=IPRIM+1
CALL NORMBF(NSHL,IPST)
DO III=0,NPRIM(NSHL)-1
WRITE(IOUQMC,'(2F20.12)') EXX(IPST+III),COEFF(IPST+III)
END DO
IPRIM=IPRIM+NPRIM(NSHL)
END DO
END DO
WRITE(IINQMC,8307)
8307 FORMAT('***** parameter rlambij ****** ',/,' 1.0')
WRITE(IOUQMC,*)
WRITE(IOUQMC,*) ' SPECIFIC NORMATION FACTORS: '
WRITE(IOUQMC,*) ' L M FACTOR FUNCTION '
DO L=1,LMAX
DO M=-L+1,L-1
XNORM(L,M)=FAC(L-ABS(M))/FAC(L+ABS(M))
c$$$ WRITE(IOUQMC,*) ' L, M, (L-|M|)!, (L+|M|)! ',L-1,M, FAC(L-ABS(M))
c$$$ $ ,FAC(L+ABS(M))
IF (M.NE.0) XNORM(L,M)=XNORM(L,M)*2.0D0
XNORM(L,M)=XNORM(L,M)*XCOEF(L,M)
WRITE(IOUQMC,8006) L-1,M,XNORM(L,M),CSTRI(L,M)
END DO
END DO
8006 FORMAT(2I3,F20.12,' ',A)
WRITE(IOUQMC,*)
C
WRITE(6,*) ' NBAS =',NBASY
WRITE(6,*) ' NOCC =',NOCC
WRITE(6,*) ' LMAX =',LMAX
WRITE(6,*) ' NSHL =',NSHL
WRITE(6,*) ' NPRIM =',(NPRIM(I),I=1,NSHL)
WRITE(6,*) ' NPRIMT =',NPRIMT
C
WRITE(IINQMC,8103) 2*NOCC,NOCC,NOCC
WRITE(IINQMC,8104)
WRITE(IINQMC,8105) NUMDET
8103 FORMAT('Number of electrons',/,I4,/
$ ,'Number of elect. alpha and beta',/,2I4)
8104 FORMAT('************************************',/
$ ,'Begin: Wave function',/
$ ,'***************************************************',/
$ ,'Simple Jast factor: enter 1; Sophis. Jast : enter 0',/,'1',
$ /,'************************************',/
$ ,'Pseudo: enter 1; No Pseudos : enter 0',/,'0',/
$ ,'***************************************************')
8105 FORMAT('Number of configurations:',/,I4,/
$ ,'******************************************',/
$ ,'STARTING READING PARAMETERS SIMPLE JASTROW',/
$ ,'******************************************',/
$ ,'Parameters of simple Jastrow function: ')
8106 FORMAT(2I3,F10.6)
8107 FORMAT('Number of parameters:',/,I4,/
$ ,'END READING PARAMETERS SIMPLE JASTROW',/
$ ,'**************************************************',/
$ ,'BEGINNING READING PARAMETERS SOPHISTICATED JASTROW',/
$ ,'**************************************************',/
$ ,'Parameters of sophisticate Jastrow function:')
C
NJASIM=4+NATOM
NJASOF=23+5*NATOM
IZERO=0
ZERO=0.D0
DO IDET=1,NUMDET
DO IPAR=1,NJASIM
WRITE(IINQMC,8106) IPAR,IZERO,ZERO
END DO
END DO
WRITE(IINQMC,8107) NJASIM
DO IDET=1,NUMDET
DO IPAR=1,NJASOF
WRITE(IINQMC,8106) IPAR,IZERO,ZERO
END DO
END DO
WRITE(IINQMC,8108) NJASOF
8108 FORMAT('Number of parameters',/,I4,/
$ ,'END READING PARAMETERS SOPHISTICATED JASTROW',/
$ ,'**********************************',/,'Theta',/,'0.',/
$ ,'Atomic basis functions with or without normalization'
$ ,' prefactor',/
$ ,'1----> with prefactor 0-----> without prefactor',/,'0')
C
C that is all for the zeroed Jastrow
C
WRITE(6,*)
WRITE(IINQMC,8109) NPRIMT
8109 FORMAT('Nombre de fonctions de base:',/,I5,/
$ ,'Center Type Case exponent')
8110 FORMAT(3I6,F20.12,' ',A16)
C
C reconstruct the angular momenta of the basis functions
C
C we assume that they are ordered with respect to angular momenta
C
WRITE(6,*) ' THE BASIS FUNCTIONS AND THEIR ANGULAR MOMENTA '
WRITE(6,*)
WRITE(6,*) ' IBAS CENTER L M '
WRITE(6,*) ' ============================'
IBAS=0
DO IAT=1,NATOM
DO ILL=1,LMAX
DO ILY=1,NSHLT(ILL,IAT)
C p functions are x y z
IF (ILL.EQ.2) THEN
IKK=1
IBAS=IBAS+1
WRITE(6,9221) IBAS,IAT,CST(ILL),IKK
IKK=-1
IBAS=IBAS+1
WRITE(6,9221) IBAS,IAT,CST(ILL),IKK
IKK=0
IBAS=IBAS+1
WRITE(6,9221) IBAS,IAT,CST(ILL),IKK
ELSE
DO IKK=-ILL+1,ILL-1
IBAS=IBAS+1
WRITE(6,9221) IBAS,IAT,CST(ILL),IKK
END DO
END IF
END DO
END DO
END DO
WRITE(6,*)
$ ' ====================================================== '
WRITE(6,*)
9221 FORMAT(2X,I6,I7,4X,A1,I6)
WRITE(6,*)
WRITE(6,*)
WRITE(6,*) 'AT THE LEVEL OF PRIMITIVES: '
WRITE(6,*)
WRITE(6,*) ' IPCOUNT IPRIM CENTR L '
$ ,'M Exp old Coeff new Coeff'
WRITE(6,*) ' ================================'
$ ,'========================================= '
C
WRITE(IOUQMC,*)
WRITE(IOUQMC,*)
$ ' IPRIM CENTER L M EXPONENT TYPE '
IPRIM=0
ISHL=0
IPCOUNT=0
IBAS=0
C..FILE 'morceau1.h'
DO IAT=1,NATOM
DO ILL=1,LMAX
DO ILY=1,NSHLT(ILL,IAT)
ISHL=ISHL+1
LVAL=IL(ISHL)
NPRIML=NPRIM(ISHL)
IPST=IPRIM
C p functions, again
IF (ILL.EQ.2) THEN
IKK=1
IPRIM=IPST
IBAS=IBAS+1
DO III=1,NPRIML
IPRIM=IPRIM+1
EXPO=EXX(IPRIM)
IPCOUNT=IPCOUNT+1
WRITE(6,9222) IPCOUNT,IPRIM,IAT,LVAL,IKK,EXPO
$ ,COEFF(IPRIM),COEFF(IPRIM)*XNORM(LVAL+1,IKK)
WRITE(IOUQMC,9223) IPCOUNT,IAT,LVAL,IKK,EXPO,
- CSTRI(LVAL+1,IKK)
WRITE(IINQMC,8110) IAT,LVAL+9,IKK,EXPO,CSTRI(LVAL+1,IKK)
CALL NORM_PRIM(EXPO,LVAL,IKK,DCOEFF)
XNBUF(IPCOUNT)=COEFF(IPRIM)*DCOEFF*XCOEF(LVAL+1,IKK)
IF (ABS(XNBUF(IPCOUNT)).LE.1.D-6) WRITE(6,*) ' DOUBTFUL :'
$ ,COEFF(IPRIM),DCOEFF,XCOEF(LVAL+1,IKK)
IWBUF(IPCOUNT)=IBAS
END DO
C.. INCLUDE 'morceau1.h'
IKK=-1
IPRIM=IPST
IBAS=IBAS+1
DO III=1,NPRIML
IPRIM=IPRIM+1
EXPO=EXX(IPRIM)
IPCOUNT=IPCOUNT+1
WRITE(6,9222) IPCOUNT,IPRIM,IAT,LVAL,IKK,EXPO
$ ,COEFF(IPRIM),COEFF(IPRIM)*XNORM(LVAL+1,IKK)
WRITE(IOUQMC,9223) IPCOUNT,IAT,LVAL,IKK,EXPO,
- CSTRI(LVAL+1,IKK)
WRITE(IINQMC,8110) IAT,LVAL+9,IKK,EXPO,CSTRI(LVAL+1,IKK)
CALL NORM_PRIM(EXPO,LVAL,IKK,DCOEFF)
XNBUF(IPCOUNT)=COEFF(IPRIM)*DCOEFF*XCOEF(LVAL+1,IKK)
IF (ABS(XNBUF(IPCOUNT)).LE.1.D-6) WRITE(6,*) ' DOUBTFUL :'
$ ,COEFF(IPRIM),DCOEFF,XCOEF(LVAL+1,IKK)
IWBUF(IPCOUNT)=IBAS
END DO
C.. INCLUDE 'morceau1.h'
IKK=0
IPRIM=IPST
IBAS=IBAS+1
DO III=1,NPRIML
IPRIM=IPRIM+1
EXPO=EXX(IPRIM)
IPCOUNT=IPCOUNT+1
WRITE(6,9222) IPCOUNT,IPRIM,IAT,LVAL,IKK,EXPO
$ ,COEFF(IPRIM),COEFF(IPRIM)*XNORM(LVAL+1,IKK)
WRITE(IOUQMC,9223) IPCOUNT,IAT,LVAL,IKK,EXPO,
- CSTRI(LVAL+1,IKK)
WRITE(IINQMC,8110) IAT,LVAL+9,IKK,EXPO,CSTRI(LVAL+1,IKK)
CALL NORM_PRIM(EXPO,LVAL,IKK,DCOEFF)
XNBUF(IPCOUNT)=COEFF(IPRIM)*DCOEFF*XCOEF(LVAL+1,IKK)
IF (ABS(XNBUF(IPCOUNT)).LE.1.D-6) WRITE(6,*) ' DOUBTFUL :'
$ ,COEFF(IPRIM),DCOEFF,XCOEF(LVAL+1,IKK)
IWBUF(IPCOUNT)=IBAS
END DO
C.. INCLUDE 'morceau1.h'
ELSE
DO IKK=-ILL+1,ILL-1
IPRIM=IPST
IBAS=IBAS+1
DO III=1,NPRIML
IPRIM=IPRIM+1
EXPO=EXX(IPRIM)
IPCOUNT=IPCOUNT+1
WRITE(6,9222) IPCOUNT,IPRIM,IAT,LVAL,IKK,EXPO
$ ,COEFF(IPRIM),COEFF(IPRIM)*XNORM(LVAL+1,IKK)
WRITE(IOUQMC,9223) IPCOUNT,IAT,LVAL,IKK,EXPO,
- CSTRI(LVAL+1,IKK)
WRITE(IINQMC,8110) IAT,LVAL+9,IKK,EXPO,CSTRI(LVAL+1,IKK)
CALL NORM_PRIM(EXPO,LVAL,IKK,DCOEFF)
XNBUF(IPCOUNT)=COEFF(IPRIM)*DCOEFF*XCOEF(LVAL+1,IKK)
IF (ABS(XNBUF(IPCOUNT)).LE.1.D-6) WRITE(6,*) ' DOUBTFUL :'
$ ,COEFF(IPRIM),DCOEFF,XCOEF(LVAL+1,IKK)
IWBUF(IPCOUNT)=IBAS
END DO
C.. INCLUDE 'morceau1.h'
END DO
END IF
C WRITE(6,*) ' A LA FIN :',IPRIM,IPST+NPRIML
IPRIM=IPST+NPRIML
END DO
END DO
END DO
9222 FORMAT(I6,I6,I7,I4,I3,F15.8,F12.8,' ->',F12.8)
9223 FORMAT(I6,I7,I4,I3,F20.12,' ',A)
9224 FORMAT(I6,I5,I6,I5,I3,' ',A)
C
WRITE(6,*)
WRITE(6,*) ' NORMALIZATION FACTORS OF THE INDIVIDUAL PRIMITIVES: '
WRITE(6,*) ' Primitive Basis multiplicator '
DO I=1,IPCOUNT
WRITE(6,'(I10,I7,F20.12)') I,IWBUF(I),XNBUF(I)
END DO
WRITE(6,*)
WRITE(6,*)
WRITE(6,*)
$ ' ===================================================='
WRITE(6,*)
WRITE(6,*) ' COUNTING M VALUES WE HAVE ',IPCOUNT,
- ' PRIMITIVE BASIS FUNCTIONS'
WRITE(6,*)
C
C we put down the list of determinants
C
INUM=0
IZERO=0
ONE=1.D0
IONE=1
WRITE(IOUQMC,8008) IONE,(IZERO,J=1,4),ONE
WRITE(6,8008) IONE,(IZERO,J=1,4),ONE
WRITE(6,*)
8008 FORMAT(' ',I8,': (',2I5,' ) --> (',2I5,' ) COEFF: ',F20.12)
C
WRITE(IOUQMC,8009) ETOT
8009 FORMAT(//,' Hartree-Fock Energy: ',/,F20.12,//)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
WRITE(6,*)
WRITE(6,*) ' THE REFERENCE '
WRITE(6,*)
$ ' ========================================================= '
WRITE(6,*) ' ORBITALS:'
CALL PFUNC(CI,NOCC,NBAS)
WRITE(6,*)
WRITE(6,*) ' EXPANDED ORBITALS:'
DO IORB=1,NOCC
CALL EXPNND(CI(1,IORB),OWBUF,IWBUF,XNBUF,IPCOUNT)
WRITE(6,9901) IORB,'alpha spin',(OWBUF(IPRIM),IPRIM=1,IPCOUNT)
END DO
DO IORB=1,NOCC
CALL EXPNND(CI(1,IORB),OWBUF,IWBUF,XNBUF,IPCOUNT)
WRITE(6,9901) IORB,'beta spin ',(OWBUF(IPRIM),IPRIM=1,IPCOUNT)
END DO
WRITE(6,*)
INUM=1
ONE=1.D0
WRITE(IOUQMC,9000) INUM,ONE
INDXW=INUM+NJASIM
WRITE(IINQMC,8111) INUM,INUM+NJASIM,IZERO,ONE
9000 FORMAT(' Determinant ',I5,' Weight and coefficients ',/,F20.12)
C
IDET=1
WRITE(6,*) ' WRITING DETERMINANT ',IDET,' ON FILE ... '
8111 FORMAT('Determinant ',I4,' Weight and coefficients ',/,
- 2I4,F20.12)
C
DO II=1,NOCC
IVALPH(II)=II
IVBETA(II)=II
END DO
C
IZERO=0
DO IORB=1,NOCC
CALL EXPNND(CI(1,IVALPH(IORB)),OWBUF,IWBUF,XNBUF,IPCOUNT)
WRITE(IOUQMC,9901) IORB,'alpha spin',(OWBUF(IPRIM)
$ ,IPRIM=1,IPCOUNT)
C//SKIPCN
cN WRITE(IINQMC,8112) IORB,'alpha spin',(OWBUF(IPRIM)
cN $ ,IPRIM=1,IPCOUNT)
cN/ENDCN
WRITE(IINQMC,8312) IORB,'alpha spin'
DO IPRIM=1,IPCOUNT
INDXW=INDXW+1
WRITE(IINQMC,'(2I6,D20.12)') INDXW,IZERO,OWBUF(IPRIM)
END DO
END DO
DO IORB=1,NOCC
CALL EXPNND(CI(1,IVBETA(IORB)),OWBUF,IWBUF,XNBUF,IPCOUNT)
WRITE(IOUQMC,9901) IORB,'beta spin',(OWBUF(IPRIM)
$ ,IPRIM=1,IPCOUNT)
C//SKIPCN
cN WRITE(IINQMC,8112) IORB,'beta spin',(OWBUF(IPRIM)
cN $ ,IPRIM=1,IPCOUNT)
cN/ENDCN
WRITE(IINQMC,8312) IORB,'beta spin'
DO IPRIM=1,IPCOUNT
INDXW=INDXW+1
WRITE(IINQMC,'(2I6,D20.12)') INDXW,IZERO,OWBUF(IPRIM)
END DO
END DO
9901 FORMAT(' Molecular orbital No ',I3,' ',A,/,3(F20.12))
8112 FORMAT(' Molecular orbital No ',I3,' ',A,/,3(F23.17))
8312 FORMAT(' Molecular orbital No ',I3,' ',A)
CLOSE(IOUQMC)
WRITE(IINQMC,8313)
8313 FORMAT('*** Observable part. iobsread=1 --> reading'
$ ,' qmcinobs npobs=nb of parameters ',/,'0 0',/
$ ,'*** End of observable part')
C
CLOSE(IINQMC)
C
RETURN
END
C
SUBROUTINE EXPNND(CI,OWBUF,IWBUF,XNBUF,IPCOUNT)
INCLUDE 'param.h'
DIMENSION CI(*),IWBUF(IPCOUNT),XNBUF(IPCOUNT),OWBUF(IPCOUNT)
C
DO I=1,IPCOUNT
OWBUF(I)=CI(IWBUF(I))*XNBUF(I)
END DO
RETURN
END
C
SUBROUTINE NORM_PRIM(EXPON,LVAL,MVAL,FACTOR)
INCLUDE 'param.h'
COMMON /POLYTA/ FAC(NLMAX+NLMAX),FACC(NLMAX)
C
C here we normalize a primitive Gaussian function
C
C sqrt[
C l-dependent:
C 2^l*(2a)^(l+3/2)/(pi^3/2 (2l-1)!!)
C m-dependent:
C [(l-|m|)!/(l+|m|)!] * (2-delta_m0)
C ]
C
PI=2.0D0*ACOS(0.D0)
XVAL=DBLE(LVAL)
PIF=FACC(LVAL+1)*(PI**1.5D0)
CD WRITE(6,*) ' FACC, PI ',FACC(LVAL+1),PI
PIF=(2.D0**XVAL)/PIF
CCC=(2.D0*EXPON)**(XVAL+1.5D0)
XMVAL=FAC(LVAL-ABS(MVAL)+1)/FAC(LVAL+ABS(MVAL)+1)
IF (MVAL.NE.0) XMVAL=XMVAL*2.D0
FACTOR=SQRT(PIF*CCC*XMVAL)
CD WRITE(6,*) ' NORM_PRIM: L,M,EXPON,XMVAL,FACTOR ',LVAL,MVAL,EXPON
CD $ ,XMVAL,FACTOR
RETURN
END
C
SUBROUTINE NORMBF(NSHLF,IPST)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /POLYTA/ FAC(NLMAX+NLMAX),FACC(NLMAX)
C the (2l+1)!! and the (l+|m|)!
C
C we normalize the functions x^l \sum_i c_i exp{-a_i * (r^2))} to unity
C this gives as factor:
C
C sqrt[(\pi^{3/2}*(2l-1)!! / 2^l) * \sum_{ij} c_i c_j / ((a_i + a_j)^(l+3/2))]
C
C in general this will be multiplied by
C
C sqrt[ (l+|m|)! / ((2-\delta_{m0})*(l-|m|)! ) ], a factor of 1 for m=0
C
C this is the prefactor for the solid harmonics
C
C
C well, we normalize here primitives to 1 and contractions modulo the
c individual norm of the primitives
C
XNRM=0.D0
C 0: s, 1: p, ...
LVAL=IL(NSHLF)
NPRIMM=NPRIM(NSHLF)
EPS=1.D-9
PI=2.0D0*ACOS(0.D0)
XVAL=DBLE(LVAL)
LVAL=LVAL+1
PIF=FACC(LVAL)*(PI**1.5D0)
PIF=(2.D0**XVAL)/PIF
CD WRITE(6,*) 'NORMBF',LVAL,IPST,NPRIM,EXX(IPST),FACC(LVAL)
C
IF (NPRIMM.EQ.1) THEN
CD CCC=(2.D0*EXX(IPST))**(XVAL+1.5D0)
CD COEFF(IPST)=SQRT(PIF*CCC)
COEFF(IPST)=1.D0
CD WRITE(6,*) ' a) LVAL,COEFF,EXP',IPST,LVAL-1,COEFF(IPST),EXX(IPST)
ELSE
CCC=0.D0
DO IP=IPST,IPST+NPRIMM-1
CIP=COEFF(IP)
DIP=SQRT((2.D0*EXX(IP))**(XVAL+1.5D0))
DO JP=IPST,IPST+NPRIMM-1
EXXX=1.D0/(EXX(IP)+EXX(JP))
CJP=COEFF(JP)
DJP=SQRT((2.D0*EXX(JP))**(XVAL+1.5D0))
CCC=CCC+CIP*CJP*DIP*DJP*(EXXX**(XVAL+1.5D0))
END DO
END DO
DO IP=IPST,IPST+NPRIMM-1
COEFF(IP)=COEFF(IP)/SQRT(CCC)
CD WRITE(6,*) ' b) LVAL,COEFF,EXP',IP,LVAL-1,COEFF(IP),EXX(IP)
END DO
END IF
RETURN
END
C
SUBROUTINE FFCAL
INCLUDE 'param.h'
COMMON /POLYTA/ FAC(NLMAX+NLMAX),FACC(NLMAX)
C DELCD
C
C the (2l-1)!! and the (l+|m|)!
C first the (2l-1)!!
C
FACC(1)=1.D0
DO L=2,NLMAX
FACC(L)=FACC(L-1)*DBLE(L+L-3)
END DO
C
FAC(1)=1.D0
FAC(2)=1.D0
DO I=3,NLMAX+NLMAX
XMLT=DBLE(I-1)
FAC(I)=XMLT*FAC(I-1)
END DO
C
CD DO L=1,NLMAX
CD WRITE(6,*) ' L, FAC(L), FACC(L) ',L-1,FAC(L),FACC(L)
CD END DO
CD DO L=1+NLMAX,NLMAX+NLMAX
CD WRITE(6,*) ' L, FAC(L) ',L-1,FAC(L)
CD END DO
C
RETURN
END
C
SUBROUTINE TRANSDA
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
C
LOGICAL LCORE,LREAD
DIMENSION TRFTMP(NBASM,NBASM)
C
C this routine is the DA routine from vind.r, not taking at all into
C account localization and/or linear scaling
C
WRITE(6,*)
WRITE(6,*) ' partial four-index transformation '
WRITE(6,*)
C do we store all in-core?
C
C count integrals
C
c$$$ WRITE(6,*) ' SIZE PARAMETERS '
c$$$ WRITE(6,*) ' NBAO = ',NBAO
c$$$ WRITE(6,*) ' NBMO = ',NBMO
C
NTRIA=NBAS*(NBAS+1)/2
C
C we know that there are less MO integrals than AO integrals, but it is
C a little complicated
C
NTRIB=NTRIA
C in AOs
NAO=NTRIA*(NTRIA+1)/2
NAOTMP=NTRIA*NTRIA
C half-transformed
NHFTR=NTRIA*NTRIB
C in MOs
NMO=NTRIB*(NTRIB+1)/2
C
C we use only the buffer NBULL, and one temporary array for the blocks
C
NMX=MAX(NAO,NAOTMP)
NMX=MAX(NMX,NHTRF)
NMX=MAX(NMO,NMX)
C
LCORE=NMX.LE.NBULL
C
IF (LCORE) THEN
C the completely-in-core transformation
C
WRITE(6,*) ' DOING THE WHOLE TRANSFORMATION IN-CORE'
WRITE(6,*)
DO I=1,NAOTMP
H0(I)=0.D0
END DO
C
IWFILE=4
CALL WOPEN(IWFILE,2)
IRD=0
NJET=0
100 CONTINUE
CALL WREAD(LREAD)
C
DO III=1,INTRD
IRD=IRD+1
HHH=HBUF(III)
C
IF (ABS(HHH).LT.THRESH) THEN
NJET=NJET+1
ELSE
I1=IBUF(III)
I2=JBUF(III)
I3=KBUF(III)
I4=LBUF(III)
C
CALL OCAN2(I1,I2,I3,I4,WEIGHT)
C
C IDBTRI gives the position in a double triangle
C
IFI34=IDBTRI(I1,I2,I3,I4,NBAS)
IFI12=IDBTRI(I3,I4,I1,I2,NBAS)
C
H0(IFI34)=HHH*WEIGHT
H0(IFI12)=HHH*WEIGHT
C
END IF
C
C close loop over integrals in one buffer
C
END DO
IF (.NOT.LREAD) GO TO 100
C
CALL WCLOS(2)
C
WRITE(6,*) ' READ ',IRD,' BIELECTRONIC INTEGRALS '
WRITE(6,*) ' REJECTED ',NJET,' INTEGRALS '
WRITE(6,*)
C
C all AO integrals are in-core, in blocks of length NTRIA
C
INDX0=0
DO IALPH=1,NBAS
DO IBETA=IALPH,NBAS
INDX1=INDX0
DO IGAMM=1,NBAS
DO IDELT=IGAMM,NBAS
INDX1=INDX1+1
HDUM=H0(INDX1)
TRFTMP(IGAMM,IDELT)=HDUM
TRFTMP(IDELT,IGAMM)=HDUM
END DO
END DO
C
C (ab|cd) -> (on|cd)
C
INDX1S=1
INDX1F=NOCC
INDX2S=1
INDX2F=NBAS
CALL TRANSL(TRFTMP,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
C
C copy the result into H0, at the double-triangle positions
C
INDX1=INDX0
DO IC=1,NOCC
DO ID=IC,NBAS
INDX1=INDX0+ITRANG(IC,ID,NBAS)
H0(INDX1)=TRFTMP(IC,ID)
END DO
END DO
INDX0=INDX0+NTRIA
END DO
END DO
C
C now the second half-transformation
C
DO IC=1,NOCC
DO ID=IC,NBAS
INDX0=ITRANG(IC,ID,NBAS)
INDX1=INDX0
DO IALPH=1,NBAS
DO IBETA=IALPH,NBAS
HDUM=H0(INDX1)
TRFTMP(IALPH,IBETA)=HDUM
TRFTMP(IBETA,IALPH)=HDUM
INDX1=INDX1+NTRIA
END DO
END DO
C
C (on|cd) -> (on|mv)
C
INDX1S=1
INDX1F=NBAS
INDX2S=NOCC+1
INDX2F=NBAS
CALL TRANSL(TRFTMP,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
C
C we could store the transformed integrals right away onto file
C
INDX1=INDX0
DO IA=1,NBAS
DO IB=IA,NBAS
C store only the mv integrals in H0
IF (IA.GT.NOCC+1) H0(INDX1)=TRFTMP(IA,IB)
INDX1=INDX1+NTRIA
END DO
END DO
C
END DO
END DO
C
C we contract: 2*(ia|jb)-(ij|ab), and output to file
C
C we write the integrals onto file 12
IUNITW=12
CALL WOPEN(IUNITW,1)
DO I=1,NOCC
DO IA=NOCC+1,NBAS
DO J=I+1,NOCC
DO IB=NOCC+1,NBAS
IF (IA.NE.IB) THEN
C integral (ia|jb)
IAJB=IDBTRI(I,IA,J,IB,NBAS)
C integral (ij|ab)
IAA=MIN(IA,IB)
IBB=MAX(IA,IB)
IJAB=IDBTRI(I,J,IAA,IBB,NBAS)
C
HHH=2.D0*H0(IAJB)-H0(IJAB)
CALL WADD(I,IA,J,IB,HHH)
END IF
C
END DO
END DO
END DO
END DO
CALL WCLOS(1)
C
C that was the direct in-core transformation
C
ELSE
C
C this is the out-of-core transformation
C
C how many couples gamma,delta fit in the buffer?
NCDMX=NBULL/NTRIA
NCORL=NTRIA/NCDMX+1
WRITE(6,*) ' NBULL = ',NBULL
WRITE(6,*)
WRITE(6,*) ' we need ',NCORL,
- ' passages of the bielectronic integrals'
WRITE(6,*) ' each passage we pick max. ',NCDMX
- ,' couples gamma/delta'
WRITE(6,*)
c$$$ WRITE(6,*) ' this makes ',NCDMX*NCORL,' index couples '
c$$$ WRITE(6,*) ' we actually need ',NTRIA,' index couples '
C
WRITE(6,*) ' we fill ',DBLE(NCDMX*NTRIA)/DBLE(NBULL)
- ,'% of the main buffer '
WRITE(6,*) ' that is ',NCDMX*NTRIA,' of the ',NBULL
- ,' places, leaving ',NBULL-NCDMX*NTRIA,' places empty '
WRITE(6,*)
IWFILE=7
IRD=0
NJET=0
C
IUNITR=4
C
WRITE(6,*) ' DA-File has a record length of',NCDMX
IUNIT1=29
OPEN(UNIT=IUNIT1,FILE='DAFILE.TMP',STATUS='UNKNOWN',
- ACCESS='DIRECT',RECL=NCDMX*8)
CLOSE(IUNIT1,STATUS='DELETE')
OPEN(UNIT=IUNIT1,FILE='DAFILE.TMP',STATUS='NEW',
- ACCESS='DIRECT',RECL=NCDMX*8)
NCDGO=NTRIA
ITRMN=1
ITRMX=NCDMX
IDREC=0
DO ICOREL=1,NCORL
CALL WOPEN(IUNITR,2)
NCD=MIN(NCDMX,NCDGO)
WRITE(6,*) ' read No ',ICOREL,' of integrals:'
$ ,NCD,' couples gamma/delta '
DO I=1,NCD*NTRIA
H0(I)=0.D0
END DO
C
101 CONTINUE
CALL WREAD(LREAD)
C
DO III=1,INTRD
IRD=IRD+1
HHH=HBUF(III)
C
IF (ABS(HHH).LT.THRESH) THEN
NJET=NJET+1
ELSE
I1=IBUF(III)
I2=JBUF(III)
I3=KBUF(III)
I4=LBUF(III)
C
CALL OCAN2(I1,I2,I3,I4,WEIGHT)
C
C which block?
IBLK1=ITRANG(I3,I4,NBAS)
IBLK2=ITRANG(I1,I2,NBAS)
C
IF (IBLK1.GE.ITRMN.AND.IBLK1.LE.ITRMX) THEN
INDX=(IBLK1-ITRMN)*NTRIA+IBLK2
H0(INDX)=HHH*WEIGHT
END IF
IF (IBLK2.GE.ITRMN.AND.IBLK2.LE.ITRMX) THEN
INDX=(IBLK2-ITRMN)*NTRIA+IBLK1
H0(INDX)=HHH*WEIGHT
END IF
C
END IF
C
C close loop over integrals in one buffer
C
END DO
IF (.NOT.LREAD) GO TO 101
CALL WCLOS(2)
C
C we passed once through the whole bielectronic file
C now we transform the blocks
C
INDX0=1
DO IBLK=1,NCD
INDX1=INDX0
DO IALPH=1,NBAS
DO IBETA=IALPH,NBAS
TRFTMP(IALPH,IBETA)=H0(INDX1)
TRFTMP(IBETA,IALPH)=H0(INDX1)
INDX1=INDX1+1
END DO
END DO
INDX1S=1
INDX1F=NOCC
INDX2S=1
INDX2F=NBAS
CALL TRANSL(TRFTMP,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
DO IA=1,NOCC
DO IB=IA,NBAS
INDX1=INDX0+ITRANG(IA,IB,NBAS)
H0(INDX1)=TRFTMP(IA,IB)
END DO
END DO
INDX0=INDX0+NTRIA
END DO
C
DO IA=1,NOCC
DO IB=IA,NBAS
INDX0=ITRANG(IA,IB,NBAS)
IDREC=IDREC+1
WRITE(IUNIT1,REC=IDREC) (H0(INDX0+NTRIA*(I-1)),I=1,NCDMX)
END DO
END DO
C
C close loop over integral read
NCDGO=NCDGO-NCDMX
ITRMN=ITRMX+1
ITRMX=ITRMX+NCDMX
END DO
C
WRITE(6,*)
WRITE(6,*) ' READ ',IRD,' BIELECTRONIC INTEGRALS '
WRITE(6,*) ' REJECTED ',NJET,' INTEGRALS '
WRITE(6,*)
WRITE(6,*) ' >>>> WROTE ',IDREC,' RECORDS OF LENGTH',NCDMX
WRITE(6,*)
CALL TIMING('TRF1')
C
WRITE(6,*) ' second transformation '
CLOSE(IUNIT1)
OPEN(UNIT=IUNIT1,FILE='DAFILE.TMP',STATUS='OLD',
- ACCESS='DIRECT',RECL=NCDMX*8)
C
C we open file No 13, file No 12 is for the final assembly
IUNITW=13
CALL WOPEN(IUNITW,1)
C
C again as many as possible in H0
C
C there are NTRIB blocks of size NTRIA in total
NTRIB=NOCC*(NBAS-NOCC)+NOCC*(NOCC+1)/2
NABMX=NBULL/NTRIA
WRITE(6,*) ' we may have ',NABMX,' triangles of size ',NTRIA
- ,' in the main buffer '
C and NCHNK buffer loads
NCHNK=NTRIB/NABMX+1
WRITE(6,*) ' we have ',NCHNK,' buffer loads '
C
C how many pieces of length NCDMX do we have to read
NPIECE=NTRIA/NCDMX
NREST=NTRIA-NPIECE*NCDMX
WRITE(6,*) ' each triangle gamma/delta is cut in '
- ,NPIECE,' pieces '
WRITE(6,*) ' there remain ',NREST,' index couples '
WRITE(6,*) ' this makes ',NCDMX*NPIECE+NREST
- ,' different index couples '
C
IAB=0
IREC0=1
INDX0=0
IAST=1
IBST=1
ETWO=0.D0
DO IA=1,NOCC
DO IB=IA,NBAS
IF (IAB.EQ.NABMX) THEN
WRITE(6,*)
WRITE(6,*) ' transforming and dumping the buffer: '
WRITE(6,*) ' indices start at a=',IAST,' and b=',IBST
WRITE(6,*) ' indices end at a=',IAFI,' and b=',IBFI
C in IAST, IBST we have the beginning of the chunk, and in IAFI, IBFI
C the end
WRITE(6,*) ' H0 is full up to position ',INDX0
INDX0=0
C transformation
DO IAA=IAST,IAFI
IF (IAA.EQ.IAST) THEN
IBEG=IBST
ELSE
IBEG=IAA
END IF
IF (IAA.EQ.IAFI) THEN
IBND=IBFI
ELSE
IBND=NBAS
END IF
DO IBB=IBEG,IBND
WRITE(6,*) ' A=',IAA,' B= ',IBB
- ,' treating elements ',INDX0+1,' to ',INDX0+NTRIA
C
C AO triangle
C
INDX1=INDX0
DO IGAMM=1,NBAS
DO IDELT=IGAMM,NBAS
INDX1=INDX1+1
TRFTMP(IGAMM,IDELT)=H0(INDX1)
TRFTMP(IDELT,IGAMM)=H0(INDX1)
END DO
END DO
C
INDX1S=1
INDX1F=NBAS
INDX2S=NOCC+1
INDX2F=NBAS
CALL TRANSL(TRFTMP,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
C
C IA runs over 1..NOCC, IB over IA..NBAS
C IC should run over IA..NBAS, ID over MAX(NOCC+1,IC),NBAS
C
DO IC=IAA,NBAS
IF (IC.EQ.IAA) THEN
IDST=IBB
ELSE
IDST=IC
END IF
IDST=MAX(IDST,NOCC+1)
DO ID=IDST,NBAS
HHH=TRFTMP(IC,ID)
CALL WADD(IAA,IBB,IC,ID,HHH)
END DO
END DO
INDX0=INDX0+NTRIA
END DO
END DO
WRITE(6,*) ' H0 has been transformed until position ',INDX0
IAST=IAFI
IBST=IBFI+1
IF (IBST.EQ.NBMO+1) THEN
IAST=IAST+1
IBST=IAST
END IF
IAB=0
INDX0=0
WRITE(6,*) ' ... finished '
END IF
C
IAB=IAB+1
WRITE(6,*) ' read triangle No',IAB,' for IA= ',IA,' and IB=',IB
WRITE(6,*) ' elements ',INDX0+1,' to ',INDX0+NTRIA
IREC1=IREC0
DO IPI=1,NPIECE
READ(IUNIT1,REC=IREC1,ERR=8110) (H0(INDX0+I),I=1,NCDMX)
INDX0=INDX0+NCDMX
IREC1=IREC1+NTRIB
END DO
READ(IUNIT1,REC=IREC1,ERR=8110) (H0(INDX0+I),I=1,NREST)
INDX0=INDX0+NREST
IREC0=IREC0+1
IBFI=IB
IAFI=IA
END DO
END DO
C
CLOSE(IUNIT1,STATUS='DELETE')
WRITE(6,*)
WRITE(6,*) ' deleting DA-File '
WRITE(6,*)
WRITE(6,*) ' last transformation: ',INDX0/NTRIA,' triangles '
WRITE(6,*) ' indices start at a= ',IAST,' and b= ',IBST
WRITE(6,*) ' indices end at a= ',NBMO,' and b= ',NBMO
WRITE(6,*) ' H0 is full up to position ',INDX0
C
C the last chunk to be transformed
C
INDX0=0
DO IAA=IAST,NOCC
IF (IAA.EQ.IAST) THEN
IBEG=IBST
ELSE
IBEG=IAA
END IF
DO IBB=IBEG,NBAS
INDX1=INDX0
DO IGAMM=1,NBAS
DO IDELT=IGAMM,NBAS
INDX1=INDX1+1
TRFTMP(IGAMM,IDELT)=H0(INDX1)
TRFTMP(IDELT,IGAMM)=H0(INDX1)
END DO
END DO
INDX1S=1
INDX1F=NBAS
INDX2S=NOCC+1
INDX2F=NBAS
CALL TRANSL(TRFTMP,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
DO IC=IAA,NBAS
IF (IC.EQ.IAA) THEN
IDST=IBB
ELSE
IDST=IC
END IF
IDST=MAX(IDST,NOCC+1)
DO ID=IDST,NBAS
HHH=TRFTMP(IC,ID)
CALL WADD(IAA,IBB,IC,ID,HHH)
END DO
END DO
INDX0=INDX0+NTRIA
END DO
END DO
WRITE(6,*) ' H0 has been transformed until position ',INDX0
WRITE(6,*) ' ... finished '
WRITE(6,*)
C
CALL WCLOS(1)
CALL TIMING('TRF2')
C
C now we have to contract: reread integrals (on|mv) for common ov and
C all pairs nm with n,m between o and v
C
C how many integrals do we have?
C
NTRIO=NOCC*(NOCC+1)/2
NTRIV=(NBAS-NOCC)*(NBAS-NOCC+1)/2
NOV=NOCC*(NBAS-NOCC)
NOVOV=NOV*(NOV+1)/2
NOOVV=NTRIO*NTRIV
NUMI=NOOVV+NOVOV
C
WRITE(6,*) ' we wrote ',IMOCNT,' integrals onto file out of '
$ ,NUMI,' possible integrals '
WRITE(6,*) ' now we have to construct the ',NOVOV
$ ,' integrals 2(ov|ov)-(oo|vv)'
C
C
C how many passes do we have to execute?
NPASS=NOVOV/NBULL+1
C
C we contract every integral (ov|ov) and we look for the necessary
C integrals during output
IOFFS=0
C this is the very first element
ISTRT=1
JSTRT=1
IAST=NOCC+1
IBST=IAST
C
DO IPASS=1,NPASS
DO III=1,NBULL
H0(III)=0.D0
END DO
C
CALL WOPEN(IUNITW,2)
OPEN(UNIT=12,FILE='FILE12',FORM='UNFORMATTED',STATUS='UNKNOWN')
INDXW=0
C
102 CONTINUE
CALL WREAD(LREAD)
DO III=1,INTRD
C
C we have to select and store in H0
C
C (ij|ab) -> -(ia|jb), -(ib|ja) if i.ne.j
C (ia|jb) multiply by two
C
I=IBUF(III)
J=JBUF(III)
K=KBUF(III)
L=LBUF(III)
IF (J.LE.NOCC) THEN
C
C the integral (ij|ab) goes to position (ia|jb) and to (ib|ja)
C for i.ne.j
C
IF (I.NE.J) THEN
IF (K.NE.L) THEN
INDX1=IPOSQU(I,K,J,L,NOCC,NBAS)-IOFFS
INDX2=IPOSQU(I,L,J,K,NOCC,NBAS)-IOFFS
IF (INDX1.GT.0.AND.INDX1.LE.NBULL) H0(INDX1)=H0(INDX1)-HHH
IF (INDX2.GT.0.AND.INDX2.LE.NBULL) H0(INDX1)=H0(INDX1)-HHH
END IF
END IF
ELSE
IF (I.NE.K) THEN
IF (J.NE.L) THEN
INDX=IPOSQU(I,J,K,L,NOCC,NBAS)-IOFFS
IF (INDX.GT.0.AND.INDX.LE.NBULL) H0(INDX)=H0(INDX)+2.D0*HHH
END IF
END IF
END IF
C
END DO
C
IF (.NOT.LREAD) GO TO 102
C
CALL WCLOS(2)
C write the buffer to file
C we have to include here the writing of the contracted integrals, sort
C of by hand
C
C the loop should go from (ISTRT,IAST|JSTRT,IBST) to (IFIN,IAFIN|JFIN,IBFIN)
C
I=ISTRT
J=JSTRT
IA=IAST
IB=IBST
DO III=1,NBULL
C we write the integral to file
IF (I.NE.J.AND.IA.NE.IB) THEN
INDXW=INDXW+1
IF (INDXW.EQ.IWBULL+1) THEN
WRITE(12) IBUF,JBUF,KBUF,LBUF,HBUF
INDXW=1
END IF
IBUF(INDXW)=I
JBUF(INDXW)=IA
KBUF(INDXW)=J
LBUF(INDXW)=IB
HBUF(INDXW)=H0(III)
END IF
C
CALL INCSIM(I,IA,J,IB,IERR,NOCC,NBAS)
IF (IERR.EQ.1) GO TO 2002
C
END DO
ISTRT=I
JSTRT=J
IAST=IA
IBST=IB
C
IOFFS=IOFFS+NBULL
END DO
C we jump here for the last buffer, if we have counted already all integrals
2002 CONTINUE
C write the last buffer
IF (INDXW.EQ.IWBULL) THEN
WRITE(12) IBUF,JBUF,KBUF,LBUF,HBUF
INDXW=0
END IF
IBUF(IWBULL)=-1
JBUF(IWBULL)=-1
KBUF(IWBULL)=-1
LBUF(IWBULL)=INDXW
HBUF(IWBULL)=-10000.D0
WRITE(12) IBUF,JBUF,KBUF,LBUF,HBUF
CLOSE(12)
C
CALL TIMING('CTRC')
C
END IF
C
RETURN
8100 CONTINUE
WRITE(6,*) ' NO FILE FOUND '
STOP
8110 CONTINUE
WRITE(6,*) ' Error during DA read '
WRITE(6,*) ' trying to read record No ',IREC1
WRITE(6,*) ' IREC0, IREC1, INDX0',IREC0,IREC1,INDX0
STOP
END
C
SUBROUTINE ADDMO(I,J,K,L,HHH)
INCLUDE 'param.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
C
IF (ABS(HHH).GE.THRESH) THEN
CALL WADD(I,J,K,L,HHH)
ELSE
IMOREJ=IMOREJ+1
END IF
C
RETURN
END
C
SUBROUTINE OPENMO
INCLUDE 'param.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
C
IWFILE=53
CALL WOPEN(IWFILE,1)
C
IWCNT=0
IMOCNT=0
IMOREJ=0
C
RETURN
END
C
SUBROUTINE CLOSEMO
INCLUDE 'param.h'
PARAMETER (NBULL=1000000)
C WBUF corresponds to the common block in unfor_io.r
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /TWOBUF/ H0(NBULL),IWCNT,IMOCNT,IMOREJ
C.. INCLUDE 'common_two.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
C
CALL WCLOS(1)
C
WRITE(6,*)
WRITE(6,*) ' WROTE ',IMOCNT,' BIELECTRONIC INTEGRALS '
WRITE(6,*) ' SKIPPED ',IMOREJ,' BIELECTRONIC INTEGRALS '
WRITE(6,*) ' CLOSED MO FILE '
WRITE(6,*)
C
RETURN
END
C
SUBROUTINE TRANSL(A,CI,INDX1S,INDX1F,INDX2S,INDX2F,NBAS)
INCLUDE 'param.h'
DIMENSION AJKLM(NBASM,NBASM)
DIMENSION A(NBASM,NBASM),CI(NBASM,NBASM)
C
C the two-index transformation routine, with limits for the two indices
C
DO IALPH=1,NBAS
DO J=INDX2S,INDX2F
AJKLM(I,J)=0.D0
END DO
END DO
C
DO J=INDX2S,INDX2F
DO K=1,NBAS
BB=0.D0
DO L=1,NBAS
BB=BB+CI(L,J)*A(K,L)
END DO
AJKLM(K,J)=AJKLM(K,J)+BB
END DO
END DO
C
DO J=INDX2S,INDX2F
DO I=INDX1S,INDX1F
A(I,J)=0.D0
END DO
END DO
C
DO I=INDX1S,INDX1F
DO J=INDX2S,INDX2F
SSS=0.D0
DO K=1,NBAS
SSS=SSS+AJKLM(K,J)*CI(K,I)
END DO
A(I,J)=A(I,J)+SSS
END DO
END DO
C
C copy A(I,J)=A(J,I)
C
DO I=INDX1S,INDX1F
DO J=INDX2S,INDX2F
A(J,I)=A(I,J)
END DO
END DO
C
RETURN
END
C
FUNCTION ITRANG(I,J,NBAS)
C The triangle index for J.GE.I
IF (J.LT.I) THEN
WRITE(6,*) ' I= ',I,' J=',J
STOP ' ITRANG: J Debye: ',F20.12,//
$ ,' DIPOLE MOMENT in A.U. : ',3F12.6)
9902 FORMAT(' DIPOLE MOMENT in Debye : ',3F12.6)
C
C now we do the same atom for atom as in a Mulliken analysis
C
WRITE(6,*) ' Decomposition of the total dipole '
$ ,'with a Mulliken scheme '
WRITE(6,*)
DEXT=0.D0
DEYT=0.D0
DEZT=0.D0
DO IAT=1,NATOM
WRITE(6,*) ' Atom No ',IAT
NUMAT=NZ(IAT)
DNX=NUMAT*POS(1,IAT)
DNY=NUMAT*POS(2,IAT)
DNZ=NUMAT*POS(3,IAT)
DEX=0.D0
DEY=0.D0
DEZ=0.D0
WRITE(6,*) ' nuclear',' electronic '
$ ,' total '
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
DEX=DEX+P(IALPH,IBETA)*DIPOL(IALPH,IBETA,1)
DEY=DEY+P(IALPH,IBETA)*DIPOL(IALPH,IBETA,2)
DEZ=DEZ+P(IALPH,IBETA)*DIPOL(IALPH,IBETA,3)
END DO
END DO
WRITE(6,9237) DNX,DEX,DNX+DEX,DNY,DEY,DNY+DEY,DNZ,DEZ,DNZ+DEZ
9237 FORMAT(' X ',3F15.8,/,' Y ',3F15.8,/,' Z ',3F15.8,/)
DEXT=DEXT+DEX
DEYT=DEYT+DEY
DEZT=DEZT+DEZ
END DO
RETURN
END
SUBROUTINE QPOLC
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,XFAC,ORIGIN(3),ANTOAU,XLAMB,XFMIX,LECONV,LTST
$ ,LCAN,LPIPM,LMIX,LREMO,LWFN,LCOREH,LMOLD,NITMAX,LAMBDA,LBOYS
$ ,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D,LICSCF
LOGICAL LECONV,LTST,LCAN,LPIPM,LMIX,LAMBDA,LREMO,LWFN,LCOREH,LMOLD
$ ,LBOYS,LQMC,LCICOM,LEXTER,LMULTP,LODA,LDIPOL,LREORD,L6D
$ ,LICSCF
C.. INCLUDE 'flow.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C the common things like NFRZ or NREMO we take from the Pipek-Mezey
C localization
C.. INCLUDE 'common_boys.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
COMMON /MOLPRO/ IPERM(NBASM),ISIG(NBASM),IATSRT(NATMX)
- ,ISHST(NATMX),IBFST(NATMX),CMOLP4,CMOLP7,CMOLP9,CMOLP
CHARACTER*50 CMOLP4,CMOLP7,CMOLP9,CMOLP
C.. INCLUDE 'common_bas.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION QPOLE(3,3),QN(3,3),POLINT(NBASM,NBASM),SHFTQQ(3,3)
dimension evec(3,3),eval(3),fv1(3),fv2(3),quadax(3,3,3)
DIMENSION QPMULL(6,NATMX)
C
C the quadrupole moment from the nuclei
C
DO I=1,3
DO J=1,3
QN(I,J)=0.D0
END DO
END DO
DO IAT=1,NATOM
RR=POS(1,IAT)*POS(1,IAT)+POS(2,IAT)*POS(2,IAT)+POS(3,IAT)*POS(3
$ ,IAT)
XNUMAT=0.5D0*DBLE(NZ(IAT))
DO I=1,3
QQTMP=3.D0*XNUMAT*POS(I,IAT)
DO J=1,3
QN(I,J)=QN(I,J)+QQTMP*POS(J,IAT)
END DO
QN(I,I)=QN(I,I)-XNUMAT*RR
END DO
END DO
C
C the electronic contribution to the quadrupole moment
C get the quadrupole matrix elements
C
WRITE(6,*)
WRITE(6,*) 'Calculating the quadrupole moment '
WRITE(6,*)
WRITE(6,*)
WRITE(6,9277)(ORIGIN(I),I=1,3)
9277 FORMAT(' origin: ',3F15.8)
WRITE(6,*)
RR=ORIGIN(1)*ORIGIN(1)+ORIGIN(2)*ORIGIN(2)+ORIGIN(3)*ORIGIN(3)
DO I=1,3
SHFTQQ(I,I)=0.5D0*(3.D0*ORIGIN(I)*ORIGIN(I)-RR)
DO J=1,I-1
SHFTQQ(I,J)=1.5D0*ORIGIN(I)*ORIGIN(J)
SHFTQQ(J,I)=1.5D0*ORIGIN(I)*ORIGIN(J)
END DO
END DO
WRITE(6,*) 'nuclear contribution : '
WRITE(6,9001) ((QN(I,J),I=1,3),J=1,3)
9001 FORMAT(3(15X,3F16.7,/))
WRITE(6,*)
C
C (re)generate the density matrix
DO IALPH=1,NBAS
DO IBETA=1,NBAS
P(IALPH,IBETA)=0.D0
END DO
END DO
DO IALPH=1,NBAS
DO I=1,NOCC
CCC=2.D0*CI(IALPH,I)
DO IBETA=1,NBAS
P(IALPH,IBETA)=P(IALPH,IBETA)+CCC*CI(IBETA,I)
END DO
END DO
END DO
C XX
C
C we have to shift the integrals through
C (x-x0)(x-x0)= x x - 2 x x0 + x0 x0
C (y-y0)(y-y0)= y y - 2 y y0 + y0 y0
C (x-x0)(x-x0)= z z - 2 z z0 + z0 z0
C qpole dipole overlap
C
C integrals are 1/2(3xx-rr), 3/2 xy ...
C
CALL RDMAT('OVERLAP',SAO)
CALL RDMAT('DIPOL_X',DIPOL(1,1,1))
CALL RDMAT('DIPOL_Y',DIPOL(1,1,2))
CALL RDMAT('DIPOL_Z',DIPOL(1,1,3))
CALL RDMAT('QUAD_XX',POLINT)
SHFTXX=SHFTQQ(1,1)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- - 2.D0*ORIGIN(1)*DIPOL(I,J,1)
- + ORIGIN(2)*DIPOL(I,J,2)
- + ORIGIN(3)*DIPOL(I,J,3)
- + SHFTXX*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
C
C individual contributions for each atom
C..FILE 'morceau_qq.h'
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(1,IAT)=QQQ
END DO
C
QPOLE(1,1)=QRES
C XY
SHFTXY=SHFTQQ(1,2)
CALL RDMAT('QUAD_XY',POLINT)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- - 1.5D0*(ORIGIN(1)*DIPOL(I,J,2)+ORIGIN(2)*DIPOL(I,J,1))
- + SHFTXY*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
QPOLE(1,2)=QRES
QPOLE(2,1)=QPOLE(1,2)
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(2,IAT)=QQQ
END DO
C
C XZ
CALL RDMAT('QUAD_XZ',POLINT)
SHFTXZ=SHFTQQ(1,3)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- - 1.5D0*(ORIGIN(1)*DIPOL(I,J,3)+ORIGIN(3)*DIPOL(I,J,1))
- + SHFTXZ*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
QPOLE(1,3)=QRES
QPOLE(3,1)=QPOLE(1,3)
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(3,IAT)=QQQ
END DO
C
C YY
CALL RDMAT('QUAD_YY',POLINT)
SHFTYY=SHFTQQ(2,2)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- + ORIGIN(1)*DIPOL(I,J,1)
- - 2.D0*ORIGIN(2)*DIPOL(I,J,2)
- + ORIGIN(3)*DIPOL(I,J,3)
- + SHFTYY*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
QPOLE(2,2)=QRES
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(4,IAT)=QQQ
END DO
C
C YZ
CALL RDMAT('QUAD_YZ',POLINT)
SHFTYZ=SHFTQQ(2,3)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- - 1.5D0*(ORIGIN(2)*DIPOL(I,J,3)+ORIGIN(3)*DIPOL(I,J,2))
- + SHFTYZ*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
QPOLE(2,3)=QRES
QPOLE(3,2)=QPOLE(2,3)
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(5,IAT)=QQQ
END DO
C
C ZZ
CALL RDMAT('QUAD_ZZ',POLINT)
SHFTZZ=SHFTQQ(3,3)
DO I=1,NBAS
DO J=1,NBAS
POLINT(I,J)=POLINT(I,J)
- + ORIGIN(1)*DIPOL(I,J,1)
- + ORIGIN(2)*DIPOL(I,J,2)
- - 2.D0*ORIGIN(3)*DIPOL(I,J,3)
- + SHFTZZ*SAO(I,J)
END DO
END DO
CALL DTRACE(P,POLINT,NBAS,QRES)
QPOLE(3,3)=QRES
DO IAT=1,NATOM
QQQ=0.D0
DO IALPH=IBFST(IAT),NFIN(IAT)
DO IBETA=1,NBAS
QQQ=QQQ+POLINT(IALPH,IBETA)*P(IALPH,IBETA)
END DO
END DO
C.. INCLUDE 'morceau_qq.h'
QPMULL(6,IAT)=QQQ
END DO
C
WRITE(6,*) 'electronic contribution : '
WRITE(6,9001) ((QPOLE(I,J),I=1,3),J=1,3)
WRITE(6,*)
DO I=1,3
DO J=1,3
QPOLE(I,J)=QN(I,J)-QPOLE(I,J)
END DO
END DO
C
WRITE(6,*) 'antoau = ',antoau
C
CCM = 299792458.0D0
ECHARGE = 1.602 176 462D-19
XTANG = 0.529 177 2083 D0
DEBYE = ECHARGE*XTANG*CCM*1.D11
WRITE(6,*) ' Debye of the Dalton program ',DEBYE
AUDEBYE=2.54322586D0
AUDEBYE=3.D0*1.602*antoau*antoau
C
WRITE(6,9901) AUDEBYE,((QPOLE(I,J),I=1,3),J=1,3)
WRITE(6,9902) ((QPOLE(I,J)*AUDEBYE,I=1,3),J=1,3)
WRITE(6,*)
9901 FORMAT(' conversion factor a.u. -> Debye*Angstrom: ',F20.12,//
$ ,' QUADRUPOLE MOMENT in A.U. : ',/,3(15X,3F16.8,/))
9902 FORMAT(' QUADRUPOLE MOMENT in Debye*Angstrom : ',/
$ ,3(15X,3F16.8,/))
C
C decomposition of Claverie
C
write(6,*)
write(6,*) ' Decomposition into axial quadrupoles '
write(6,*) ' - - - - - - - - - - - - - - - - - - - - - - - '
write(6,*)
CALL RS(3,3,QPOLE,EVAL,1,EVEC,FV1,FV2,IERR)
WRITE(6,9003) eval(1)
9003 FORMAT(' first eigenvalue :',F20.12)
9004 FORMAT(' second eigenvalue :',F20.12)
9005 FORMAT(' third eigenvalue :',F20.12)
WRITE(6,9002) (evec(i,1),i=1,3)
9002 FORMAT(' eigenvector ',3F16.8)
WRITE(6,*)
WRITE(6,9004) eval(2)
WRITE(6,9002) (evec(i,2),i=1,3)
WRITE(6,*)
WRITE(6,9005) eval(3)
WRITE(6,9002) (evec(i,3),i=1,3)
WRITE(6,*)
eval(3)=-eval(1)-eval(2)
evec(1,3)=evec(2,1)*evec(3,2)-evec(3,1)*evec(2,2)
evec(2,3)=evec(3,1)*evec(1,2)-evec(1,1)*evec(3,2)
evec(3,3)=evec(1,1)*evec(2,2)-evec(2,1)*evec(1,2)
write(6,*) ' reconstruction of third eigenvalue/vector :'
WRITE(6,9005) eval(3)
WRITE(6,9002) (evec(i,3),i=1,3)
WRITE(6,*)
do iq=1,3
evv=eval(iq)
do I=1,3
do j=1,3
quadax(i,j,iq)=evv*evec(i,iq)*evec(j,iq)
end do
end do
end do
write(6,*) ' first axial quadrupole : '
write(6,'(3(5X,3F16.8,/))') ((quadax(i,j,1),j=1,3),i=1,3)
WRITE(6,*)
write(6,*) ' second axial quadrupole : '
write(6,'(3(5X,3F16.8,/))') ((quadax(i,j,2),j=1,3),i=1,3)
WRITE(6,*)
write(6,*)
$ ' reassembling the axial quadrupoles to the original pole '
write(6,'(3(5X,3F16.8,/))') ((quadax(i,j,1)+quadax(i,j,2)
$ +quadax(i,j,3),j=1,3),i=1,3)
C
C Mulliken multipoles
C
WRITE(6,*) ' Decomposition of the total quadrupole '
$ ,'with a Mulliken scheme '
WRITE(6,*)
QQEXXT=0.D0
QQEXYT=0.D0
QQEXZT=0.D0
QQEYYT=0.D0
QQEYZT=0.D0
QQEZZT=0.D0
DO IAT=1,NATOM
WRITE(6,*) ' Atom No ',IAT
NUMAT=NZ(IAT)
QQNXX=NUMAT*POS(1,IAT)*POS(1,IAT)
QQNXY=NUMAT*POS(1,IAT)*POS(2,IAT)
QQNXZ=NUMAT*POS(1,IAT)*POS(3,IAT)
QQNYY=NUMAT*POS(2,IAT)*POS(2,IAT)
QQNYZ=NUMAT*POS(2,IAT)*POS(3,IAT)
QQNZZ=NUMAT*POS(3,IAT)*POS(3,IAT)
RR=QQNXX+QQNYY+QQNZZ
QQNXX= 0.5D0*(3.D0*QQNXX-RR)
QQNXY= QQNXY*1.5D0
QQNXZ= QQNXZ*1.5D0
QQNYY= 0.5D0*(3.D0*QQNYY-RR)
QQNYZ= QQNYZ*1.5D0
QQEXX=QPMULL(1,IAT)
QQEXY=QPMULL(2,IAT)
QQEXZ=QPMULL(3,IAT)
QQEYY=QPMULL(4,IAT)*1.5D0
QQEYZ=QPMULL(5,IAT)*1.5D0
QQEZZ=QPMULL(6,IAT)*1.5D0
RR=QQEXX+QQEYY+QQEZZ
QQEXX=0.5D0*(3.D0*QQEXX-RR)
QQEYY=0.5D0*(3.D0*QQEYY-RR)
WRITE(6,*) ' nuclear',' electronic '
$ ,' total '
WRITE(6,9239) QQNXX,QQEXX,QQNXX+QQEXX,QQNYY,QQEYY,QQNYY+QQEYY
$ ,QQNXY,QQEXY,QQNXY+QQEXY,QQNXZ,QQEXZ,QQNXZ+QQEXZ,QQNYZ,QQEYZ
$ ,QQNYZ+QQEYZ
9239 FORMAT(' XX ',3F15.8,/,' YY ',3F15.8,/
- ,' XY ',3F15.8,/,' XZ ',3F15.8,/,' YZ ',3F15.8,/)
DEXT=DEXT+DEX
DEYT=DEYT+DEY
DEZT=DEZT+DEZ
END DO
RETURN
END
C
SUBROUTINE DTRACE(A,B,N,Q)
INCLUDE 'param.h'
DIMENSION A(NBASM,NBASM),B(NBASM,NBASM)
C
C double trace of 2 symmetric matrices
C Q=sum_ij A_ij * B_ij = sum_i A_ii * B_ii + 2 * sum_i> molpro.o
$ut')
C folded 1 (fixf $Revision: 1.3 $)
WRITE(6,*)
WRITE(6,*) ' Molpro terminated '
WRITE(6,*)
C
C read the Fock matrix
C
IUNITP=23
OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
- FILE='gorch.mat')
READ(IUNITP,*)
C we have to pay attention to the d functions, again!!!
READ(IUNITP,*) ((F(IPERM(IALPH),IPERM(IBETA))
- ,IALPH=1,NBAS),IBETA=1,NBAS)
CLOSE(IUNITP)
C and to the sign of f0, f3
DO IALPH=1,NBAS
XMULTA=ISIG(IALPH)
DO IBETA=1,NBAS
XMULTB=ISIG(IBETA)
F(IALPH,IBETA)=XMULTA*XMULTB*F(IALPH,IBETA)
END DO
END DO
OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
- FILE='dftfun.tmp')
READ(IUNITP,*)
READ(IUNITP,*)
C read the DFT energy
READ(IUNITP,*) COUENER,EXENER,DFTENER
CLOSE(IUNITP)
C
C we cannot calculate ETWO since we do not know the functional
ETWO=COUENER+XFAC*EXENER+DFTENER
C
ELSE
C
C the bielectronic part
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
XI(IALPH,IBETA)=0.D0
END DO
END DO
IUNITR=4
CALL WOPEN(IUNITR,2)
100 CONTINUE
CALL WREAD(LREAD)
IF (LAMBDA) THEN
DO III=1,INTRD
HBUF(III)=HBUF(III)*XLAMB
END DO
END IF
C//SKIPCI
cI XNU=2.D0
cI/ENDCI
DO III=1,INTRD
I=IBUF(III)
J=JBUF(III)
K=KBUF(III)
L=LBUF(III)
HHH=HBUF(III)
C
C//SKIPCI
cI CALL ACCUG(I,J,K,L,HHH,XNU,P,XI)
cI/ENDCI
C
XI(K,L)=XI(K,L)+4.D0*P(I,J)*HHH
XI(I,J)=XI(I,J)+4.D0*P(K,L)*HHH
XI(I,K)=XI(I,K)-P(J,L)*HHH
XI(J,K)=XI(J,K)-P(I,L)*HHH
XI(I,L)=XI(I,L)-P(J,K)*HHH
XI(J,L)=XI(J,L)-P(I,K)*HHH
C
END DO
IF (.NOT.LREAD) GO TO 100
C
CALL WCLOS(2)
C
C double the matrix ... for the non-diagonal elements
DO IALPH=1,NBAS
DO IBETA=IALPH+1,NBAS
XDUM=(XI(IALPH,IBETA)+XI(IBETA,IALPH))*0.5D0
XI(IALPH,IBETA)=XDUM
XI(IBETA,IALPH)=XDUM
END DO
END DO
C
C add the simple bielectronic part to the Fock matrix
DO IALPH=1,NBAS
DO IBETA=1,NBAS
F(IALPH,IBETA)=F(IALPH,IBETA)+XI(IALPH,IBETA)
ETWO=ETWO+P(IALPH,IBETA)*XI(IALPH,IBETA)
END DO
END DO
END IF
ETWO=0.5D0*ETWO
C
c$$$C write the Fock matrix in MOLPRO format onto file
c$$$C
c$$$ IUNITP=23
c$$$ OPEN(UNIT=IUNITP,FORM='FORMATTED',STATUS='UNKNOWN',
c$$$ - FILE='FOCKM.TMP')
c$$$ WRITE(IUNITP,9235)
c$$$ 9235 FORMAT('# Fock matrix constructed by ORTHO')
c$$$ DO IALPH=1,NBAS
c$$$ WRITE(IUNITP,'(5(F15.8,1H,))') (F(IPERM(IALPH),IPERM(IBETA))
c$$$ $ ,IBETA=1,NBAS)
c$$$ END DO
c$$$ CLOSE(IUNITP)
c$$$
C
C here we arrive after the construction of the Fock matrix
C
ETOT=EONE+ETWO+EN
EELEC=EONE+ETWO
WRITE(6,*) ' electronic energy : ',EELEC
C
C save the Fock matrix for the program vind
OPEN(FILE='FOCK',FORM='FORMATTED',STATUS='UNKNOWN',UNIT=15)
WRITE(15,*) NBAS
WRITE(15,'(4E20.11)') ((F(IALPH,IBETA),IALPH=1,NBAS),IBETA=1
$ ,NBAS)
WRITE(15,*)
WRITE(15,'(4E20.11)') ETOT,EN,EONE,ETWO
CLOSE(15)
C
C the total energy
C
IF (LMULTP) THEN
WRITE(6,9902) EKIN,EHAM,EMULP,EKIN+EHAM+EMULP,ETWO,ETOT-EN
9902 FORMAT(' THE TOTAL ENERGY: ',/
$ ,' kinetic ',F20.12,/
$ ,' el-nucl ',F20.12,/
$ ,' el-multipoles ',F20.12,/
$ ,' (one-el ',F20.12,')',/
$ ,' el-el ',F20.12,/
$ ,' =====================================================',/
$ ,' total ',F20.12,' (without nuclear energy!)' ,/)
ELSE
WRITE(6,9901) EN,EKIN,EHAM,EKIN+EHAM,ETWO,ETOT,ETOT/EKIN
9901 FORMAT(' THE TOTAL ENERGY: ',/
$ ,' nuclear ',F20.12,/
$ ,' kinetic ',F20.12,/
$ ,' el-nucl ',F20.12,/
$ ,' (one-el ',F20.12,')',/
$ ,' el-el ',F20.12,/
$ ,' =====================================================',/
$ ,' total ',F20.12,' (virial theorem: ',F11.8,')',/)
END IF
RETURN
END
SUBROUTINE ICSCF
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM),NFIN(NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),XKIN(NBASM,NBASM),F(NBASM,NBASM)
$ ,HAM(NBASM,NBASM),P(NBASM,NBASM),ORBEN(NBASM),XMULP(NBASM
$ ,NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCC(NBASM),IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /ENERGY/ EN,EKIN,EONE,ETWO,ECOUL,EEXCH,ETOT
C.. INCLUDE 'common_energy.h'
DIMENSION ALPHA(NBASM,NBASM)
C
C first step: we have to get canonical orbitals
C the Fock matrix should be in molecular orbitals, canonical or local
C
C
WRITE(6,*)
WRITE(6,*) ' Internally consistent SCF -- '
$ ,'E.R.Davidson JCP 57 (1972) 1999'
WRITE(6,*)
IF (NOCC.EQ.1) THEN
WRITE(6,*)
WRITE(6,*) ' this does not work for 2 electrons, sorry '
WRITE(6,*)
RETURN
END IF
C read converged Fock matrix again
OPEN(UNIT=16,FILE='FOCK',STATUS='OLD',FORM='FORMATTED')
READ(16,*) IDUM
READ(16,*) ((F(IALPH,IBETA),IALPH= 1,NBAS),IBETA=1,NBAS)
CLOSE(16)
CALL TRANS(F,CI,NBAS,NBAS)
CALL GENCAN(F,CI,ORBEN,NBAS)
C
SUM=0.D0
DO I=1,NOCC
SUM=SUM+ORBEN(I)
END DO
C
CALL RDMAT('HAMILTO',HAM)
C transform HAM onto the canonical orbitals
CALL TRANS(HAM,CI,NBAS,NBAS)
C
C astonishingly the 2-index transformation does not give
C the same one-electron energy as the double trace of H with P
C
C so we introduce a scaling factor for correcting the error
C this is quite strange, but expolains deviations of E-4
C
EONEH2=0.D0
DO I=1,NOCC
EONEH2=EONEH2+HAM(I,I)
END DO
EONEH2=2.D0*EONEH2
WRITE(6,*) ' the sum of the orbital energies: ',2.D0*SUM
WRITE(6,*) ' one-electron error ',EONEH2-EONE
SCALE = EONEH2/EONE
C renormalize the matrix
WRITE(6,*) ' scaling factor for the H matrix ',SCALE
DO I=1,NBAS
DO J=1,NBAS
HAM(I,J)=HAM(I,J)/SCALE
END DO
END DO
EHF=0.D0
EONEH2=0.D0
DO I=1,NOCC
EHF=EHF+HAM(I,I)+ORBEN(I)
EONEH2=EONEH2+HAM(I,I)
END DO
EONEH2=2.D0*EONEH2
WRITE(6,*) ' recalculated one-electron energy ',EONEH2
WRITE(6,*) ' recalculated Hartree-Fock energy ',EHF+EN
WRITE(6,*) ' difference to previous HF energy ',EHF+EN-ETOT
WRITE(6,*) ' electronic energy ',EHF
WRITE(6,*)
NELEC=NOCC*2
WRITE(6,*) ' we have ',NELEC,' electrons '
A=0.D0
DO I=1,NOCC
A=A+ORBEN(I)-HAM(I,I)
END DO
A=A/DBLE(NELEC-1)
WRITE(6,*) ' the term A : ',A
DO I=1,NBAS
DO J=1,NBAS
ALPHA(I,J)=-HAM(I,J)
END DO
ALPHA(I,I)=ALPHA(I,I)+ORBEN(I)-A
END DO
C
C we have to rescale ALPHA
DO I=1,NOCC
DO J=1,NOCC
ALPHA(I,J)=ALPHA(I,J)/DBLE(NELEC-2)
END DO
DO J=NOCC+1,NBAS
ALPHA(I,J)=0.D0
ALPHA(J,I)=0.D0
END DO
END DO
DO I=NOCC+1,NBAS
DO J=NOCC+1,NBAS
ALPHA(I,J)=-ALPHA(I,J)/DBLE(NELEC)
END DO
END DO
TRALPH=0.5D0*A
c$$$ TRALPH=0.D0
c$$$ DO I=1,NOCC
c$$$ TRALPH=TRALPH+ALPHA(I,I)
c$$$ END DO
WRITE(6,*) ' the term tr ALPHA * RHO : ',TRALPH,0.5D0*A
C
C construct the new Fock matrix
C
DO I=1,NBAS
DO J=1,NBAS
F(I,J)=-ALPHA(I,J)
END DO
END DO
DO I=1,NBAS
F(I,I)=F(I,I)+ORBEN(I)-2.D0*TRALPH+A/NOCC
END DO
C
C and we diagonalize the new Fock matrix
C
CALL GENCAN(F,CI,ORBEN,NBAS)
C
SUM=0
DO I=1,NOCC
SUM=SUM+ORBEN(I)
WRITE(6,9901) I,ORBEN(I)
9901 FORMAT(' ICSCF Orben ',I4,F20.12)
END DO
SUM=SUM*2.D0
WRITE(6,*) ' difference to the HF electronic energy : '
$ ,SUM-EHF
WRITE(6,*)
C we should save as well the orbitals somewhere
CALL WVEC(6)
C
C the proof ?
C
C READ the Fock matrix
C OPEN(UNIT=16,FILE='FOCK',FORM='FORMATTED',STATUS='OLD')
C READ(16,*)
C READ(16,*) ((F(I,J),I=1,NBAS),J=1,NBAS)
C CLOSE(16)
C CALL TRANS(F,CI,NBAS,NBAS)
C CALL GENCAN(F,CI,ORBEN,NBAS)
C CALL RDMAT('HAMILTO',HAM)
C transform HAM onto the canonical orbitals
C CALL TRANS(HAM,CI,NBAS,NBAS)
C EHF=0.D0
C DO I=1,NOCC
C EHF=EHF+HAM(I,I)+ORBEN(I)
C END DO
C WRITE(6,*) ' recalculated Hartree-Fock energy ',EHF+EN
C WRITE(6,*) ' electronic energy ',EHF
C WRITE(6,*)
C
RETURN
END