C -*- Fortran -*-
C DELCD DEBUG
C..DELCP PROJECTION FOR INDIVI
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 DELCI inverse SMILES ordering
C
C..FILE 'flow.h'
C..FILE 'common_syst.h'
C..FILE 'common_intu.h'
C
C..FILE 'common_boys.h'
C
C..FILE 'common_extreme.h'
C
C..FILE 'common_vect.h'
C
C..FILE 'common_mezey.h'
C
PROGRAM LPMB
C
C one single molecule
C
C a posteriori localization: Boys, Pipek-Mezey, extreme
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C P.Reinhardt, 4/2003
C statistical analysis: 6/2004
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCLUDE 'param.h'
C
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(NBASM)
C.. INCLUDE 'common_mezey.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C.. INCLUDE 'common_boys.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
C
CHARACTER*8 PNAME
CHARACTER*1 BBBB
CHARACTER*4 COTYP(5)
DIMENSION ISRTI(NBASM),CIN(NBASM,NBASM)
LOGICAL LFOUND
PNAME='L P M B '
INCLUDE 'compiler_stamp'
WRITE(6,*) '$Id: lpmb.r,v 1.34 2008/04/29 09:10:30 reinh Exp $'
X=CPTIME(3)
CALL DATING(PNAME,1)
C
DATA COTYP /'Atom','Bond','Tric','Quad','Dloc'/
EPS=1.D-14
EPSE=1.D-11
NITMAX=50
C
WRITE(6,*)
WRITE(6,*) ' L P M B - Version for molecules'
WRITE(6,*) ' a posteriori localization '
WRITE(6,*)
WRITE(6,*)
WRITE(6,*) ' P.Reinhardt, Paris 4/2003 '
WRITE(6,*)
WRITE(6,*)
WRITE(6,*)
C
C read flow options from File INPUT.LPB
C
CALL RDINP
C
C print result of options scan
C
C the a posteriori (de-)localization options
C
IF (LCAN) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE CANONICAL ORBITALS '
WRITE(6,*)
ELSE
IF (LPIPM) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE PIPEK-MEZEY LOCALIZED ORBITALS '
WRITE(6,*)
ELSE
IF (LBOYS) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE BOYS LOCALIZED ORBITALS '
WRITE(6,*)
IF (LEXTR) THEN
WRITE(6,*)
WRITE(6,*) ' WE WILL PRODUCE NONORTHOGONAL, '
$ ,'EXTREMELY LOCALIZED ORBITALS '
WRITE(6,*)
ELSE
WRITE(6,*) ' NOTHING SPECIFIED: please choose either '
WRITE(6,*) ' CANONICAL, BOYS or PIPEKMEZEY in INPUT.LPB'
WRITE(6,*)
END IF
END IF
END IF
END IF
C
IF (LFRAGM) THEN
WRITE(6,*) ' FRAGMENTS defined '
ELSE
WRITE(6,*) ' explicit definition of orbital groups '
END IF
C
IF (LSTAT) THEN
WRITE(6,*) ' Statitical analysis demanded '
ELSE
WRITE(6,*) ' no statistical analysis requested '
END IF
IF (L6D) THEN
WRITE(6,*)
$ ' Orbitals will be projected onto cartesian functions '
END IF
C
IF (LSMILE) THEN
WRITE(6,*)
C--SKIPCI
WRITE(6,*) ' reordering orbital coefficients from SMILES code to
$'
C folded 1 (fixf $Revision: 1.3 $)
WRITE(6,*) ' the Dalton ordering (py,pz,px) -> (px,py,pz)'
C//ENDCI
CI WRITE(6,*) ' reordering orbital coefficients from DALTON code to '
CI WRITE(6,*) ' the SMILES ordering (px,py,pz) -> (py,pz,px)'
WRITE(6,*)
END IF
C
C read information on system
C
IUNITR=23
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
C we read the basis correctly, once and for all
C
NBAS=1
LMAX=0
NSHL=0
IPRIM=0
DO IAT=1,NATOM
READ(IUNITR,*) NZ(IAT),NSH(IAT),(POS(J,IAT),J=1,3)
DO ISH=1,NSH(IAT)
NSHL=NSHL+1
READ(IUNITR,*) ITYPE,NPRIM(NSHL)
IL(NSHL)=ITYPE
LMAX=MAX(LMAX,ITYPE+1)
DO III=1,ITYPE+ITYPE+1
C store the atomic centers of the basis functions (needed for Pipek
c -Mezey)
IATZ(NBAS)=IAT
IF (ITYPE.EQ.1) THEN
IF (III.EQ.1) THEN
C py: 1->2
CI ISMILE(NBAS+1)=NBAS
C pz: 2->3
CI ISMILE(NBAS+2)=NBAS+1
C px: 3->1
CI ISMILE(NBAS) =NBAS+2
C just the other way round
C py: 2->1
C--SKIPCI
ISMILE(NBAS) = NBAS+1
C pz: 3->2
ISMILE(NBAS+1) = NBAS+2
C px: 1->3
ISMILE(NBAS+2) = NBAS
C//ENDCI
END IF
ELSE
ISMILE(NBAS)=NBAS
END IF
IF (ITYPE.EQ.3) THEN
IF (III.EQ.4.OR.III.EQ.7) THEN
C f0 and f+2 get a sign -1
ISMILE(NBAS)=-NBAS
END IF
END IF
WRITE(6,*) ' NBAS, ITYPE = ',NBAS,ITYPE
NBAS=NBAS+1
END DO
DO III=1,NPRIM(NSHL)
IPRIM=IPRIM+1
READ(IUNITR,*) EXX(IPRIM),COEFF(IPRIM)
END DO
END DO
END DO
NBAS=NBAS-1
CLOSE(IUNITR)
C
C
IF (LSMILE) THEN
WRITE(6,*) ' the reordering of orbitals '
DO I=1,NBAS
C--SKIPCI
WRITE(6,*) ' DALTON ',I,' <- SMILES ',ISMILE(I)
C//ENDCI
CI WRITE(6,*) ' DALTON ',I,' -> SMILES ',ISMILE(I)
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
WRITE(6,*)
$ ' for Pipek-Mezey: ATOM, ISTRT, NABAS ',IAT, ISTRT(IAT),
$ NABAS(IAT)
END DO
C
IF (NBAS.GT.NBASM) THEN
WRITE(6,*) 'NBAS, MAXIMUM: ',NBAS,NBASM
STOP 'INCREASE NBASM IN param.h AND RECOMPILE'
END IF
C
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
IF (LFRAGM) THEN
WRITE(6,*)
WRITE(6,*) ' attribution of fragments: '
WRITE(6,'(6(I4,4H -> ,I3))') (IAT,IFRAGM(IAT),IAT=1,NATOM)
WRITE(6,*)
END IF
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,9301) EN
9301 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)
9902 FORMAT(/,' THE DIAGONAL ELEMENTS OF THE AO OVERLAP MATRIX'
- ,/,50(4(I5,F10.3),/))
WRITE(6,*)
WRITE(6,*) ' READ OVERLAP MATRIX INTEGRALS '
C
C read vector
C
IUNITO=31
WRITE(6,*)
WRITE(6,*) ' READING VECTOR FROM UNIT',IUNITO,' FILE: VECTOR'
OPEN(UNIT=IUNITO,FILE='VECTOR',FORM='FORMATTED',STATUS='OLD',ERR
$ =902)
READ(IUNITO,*) ((CI(I,J),I=1,NBAS),J=1,NBAS)
READ(IUNITO,*) (IOCC(I),I=1,NBAS)
READ(IUNITO,*) (IOCCS(I),I=1,NBAS)
CLOSE(IUNITO)
C
c$$$ IF (LSMILE) THEN
c$$$ DO IVEC=1,NBAS
c$$$ DO IALPH=1,NBAS
c$$$ IIPOS=ABS(ISMILE(IALPH))
c$$$ IISIG=SIGN(1,ISMILE(IALPH))
c$$$ IF (IVEC.EQ.1) WRITE(6,*) ' IIPOS, IISIG = ',IIPOS,IISIG
c$$$ CIN(IALPH,IVEC)=CI(IIPOS,IVEC)*DBLE(IISIG)
c$$$ END DO
c$$$ END DO
c$$$ DO I=1,NBAS
c$$$ DO J=1,NBAS
c$$$ CI(I,J)=CIN(I,J)
c$$$ END DO
c$$$ END DO
c$$$ END IF
NOCC=0
DO I=1,NBAS
IF (IOCCS(I).EQ.IOCC(I)) NOCC=NOCC+1
END DO
NVIRT=NBAS-NOCC
WRITE(6,*) ' NOCC ',NOCC
WRITE(6,*)
C CALL SNORM(CI,NBAS,NBAS)
C
WRITE(6,*)
WRITE(6,*) ' checking orthonormality: '
C
C using CIN as dummy array
C
DO IBETA=1,NBAS
DO IALPH=1,NBAS
CIN(IALPH,IBETA)=SAO(IALPH,IBETA)
END DO
END DO
CALL TRANS(CIN,CI,NBAS,NBAS)
OPEN(UNIT=57,FILE='OVERLAP.MO',STATUS='UNKNOWN',
- FORM='FORMATTED')
WRITE(57,'(2I5,F20.12)') ((I,J,CIN(I,J),J=1,I),I=1,NBAS)
CLOSE(57)
C
SONMX=0.D0
SOGMX=0.D0
DO I=1,NBAS
SSY=ABS(CIN(I,I)-1.D0)
IF (SSY.GT.SONMX) SONMX=SSY
DO J=I+1,NBAS
SSX=ABS(CIN(I,J))
IF (SSX.GT.SOGMX) SOGMX=SSX
END DO
END DO
C
WRITE(6,*) ' largest difference to 1 on the diagonal: ',SONMX
WRITE(6,*) ' largest off-diagonal of S: ',SOGMX
WRITE(6,*)
C
XELEC=0.D0
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
XELEC=XELEC+SAO(IALPH,IBETA)*SUM
END DO
END DO
WRITE(6,*)
WRITE(6,*) ' total number of electrons: ',2.D0*XELEC
WRITE(6,*)
C dumping the vector in the right order on a file VECTOR.SMILE for
C transferring it to a starting vector of a Dalton calculation
C
C this is not necessary if we rest inside the Slater calculation
C
IF (LSMILE) THEN
C dump the vector somewhere
WRITE(6,*)
WRITE(6,*) ' WRITING VECTOR TO '
WRITE(6,*)
IUNITX=25
OPEN(UNIT=IUNITX,FILE='VECTOR.SMILE',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
DO J=1,NBAS
WRITE(IUNITX,'(4E20.12)') (DBLE(SIGN(1,ISMILE(K)))
$ *CI(ABS(ISMILE(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)
END IF
C
C read Fock matrix
C
IUNITF=23
OPEN(UNIT=IUNITF,FILE='FOCK',FORM='FORMATTED',STATUS='OLD',ERR
$ =904)
READ(IUNITF,*) IDUM1
WRITE(6,*) ' READING FOCK MATRIX FROM UNIT',IUNITF
WRITE(6,*)
READ(IUNITF,*) ((FAO(I,J),I=1,NBAS),J=1,NBAS)
READ(IUNITF,*) EDUM,EN,E1,E2
CLOSE(IUNITF)
go to 905
904 continue
write(6,*)
$ ' no Fock matrix found, continuing without energy analysis'
DO I=1,NBAS
DO J=1,NBAS
FAO(I,J)=0.D0
END DO
END DO
EN=0.d0
E1=0.d0
E2=0.d0
905 continue
WRITE(6,*) ' E0 E1 E2 ',EN,E1,E2
WRITE(6,*)
WRITE(6,*) ' total energy: ',EN+E1+E2
WRITE(6,*)
CALL PFUNC(CI,NBAS,NBAS)
WRITE(6,*)
CALL TIMING('READ')
WRITE(6,*)
C
C we may start to work
C
DO I=1,NBAS
DO J=1,NBAS
FMO(I,J)=FAO(I,J)
END DO
END DO
C
C F in MOs
C
CALL TRANS(FMO,CI,NBAS,NBAS)
WRITE(6,*)
WRITE(6,*) ' checking Brillouin''s theorem '
FFY=0.D0
DO I=1,NOCC
DO J=NOCC+1,NBAS
IF (ABS(FMO(I,J)).GT.FFY) FFY=ABS(FMO(I,J))
END DO
END DO
WRITE(6,*) ' largest occ-virt. Fock-matrix element: ',FFY
WRITE(6,*)
C
WRITE(6,*)
C
C write out information on grouping
C
DO I=1,NBAS
ORBEN(I)=FMO(I,I)
END DO
C sort orbitals after their energy
CALL BUBBLE(NBAS,NBAS,ORBEN,FMO,CI)
C
C calulate spread of the orbitals
C
CALL SPREAD(CI)
C
C if fragments are defined we have to attribute the orbitals now
C
IF (LFRAGM) THEN
C
C how many fragments were defined?
C
IFMX=0
DO IAT=1,NATOM
IF (IFRAGM(IAT).GT.IFMX) IFMX=IFRAGM(IAT)
END DO
NGRP=IFMX*2
WRITE(6,*) ' No of fragments: ',IFMX
DO IORB=1,NBAS
IF (IOTYP(IORB).EQ.1) THEN
IF (IOCC(IORB).NE.0) THEN
NUMGRP(IORB)=IFRAGM(IATYP(1,IORB))
ELSE
NUMGRP(IORB)=IFMX+IFRAGM(IATYP(1,IORB))
END IF
ELSE
C are multi-center orbitals on the same fragment?
LFOUND=.TRUE.
DO I1=1,IOTYP(IORB)-1
IFT1=IFRAGM(IATYP(I1,IORB))
DO I2=I1+1,IOTYP(IORB)
IFT2=IFRAGM(IATYP(I2,IORB))
IF (IFT2.NE.IFT1) LFOUND=.FALSE.
END DO
END DO
IF (.NOT.LFOUND) THEN
NGRP=NGRP+1
NUMGRP(IORB)=NGRP
ELSE
IF (IOCC(IORB).NE.0) THEN
NUMGRP(IORB)=IFRAGM(IATYP(1,IORB))
ELSE
NUMGRP(IORB)=IFMX+IFRAGM(IATYP(1,IORB))
END IF
END IF
END IF
END DO
END IF
WRITE(6,9400)
9400 FORMAT(//,7X,'No',11X,'orb.-en.',14X,'occ. group type')
DO I=1,NBAS
IF (IOTYP(I).EQ.5) THEN
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
ELSE
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
$ ,(IATYP(J,I),J=1,IOTYP(I))
END IF
END DO
9401 FORMAT(5x,I5,5x,F16.8,10X,2I5,5X,A4,5I3)
C
C we reorder for having the groups in ascending order, and within each
C group, we have the orbitals sorted by energy
C
WRITE(6,*)
CALL TIMING('PREP')
WRITE(6,*)
IF (NGRP.NE.0.OR.LEXTR) THEN
C
C will we mix occupied and virtuals?
C
INDX=0
DO IGRP=1,NGRP
CD WRITE(6,*) ' GROUP No: ',IGRP
ISTGR(IGRP)=INDX+1
IF (IGRP.GE.2) IFIGR(IGRP-1)=INDX
LFOUND=.FALSE.
DO IORB=1,NBAS
IF (NUMGRP(IORB).EQ.IGRP) THEN
IF (LFOUND) THEN
IF (IOCC(IORB).NE.IOCCC) WRITE(6,*) ' WARNING: mixing ',
$ ' occupied and virtuals !!'
ELSE
LFOUND=.TRUE.
IOCCC=IOCC(IORB)
END IF
INDX=INDX+1
ISRTI(INDX)=IORB
CD WRITE(6,*) ' Orbital No ',IORB,' goes to place ',INDX
END IF
END DO
END DO
IFIGR(NGRP)=INDX
C
DO IGRP=1,NGRP
IF (ISTGR(IGRP).LE.IFIGR(IGRP)) THEN
WRITE(6,*) ' Group No ',IGRP,' goes from ',ISTGR(IGRP),' to '
$ ,IFIGR(IGRP)
ELSE
WRITE(6,*) ' No member in Group No ',IGRP
END IF
END DO
C
DO IORB=1,NBAS
IF (NUMGRP(IORB).EQ.0) THEN
INDX=INDX+1
ISRTI(INDX)=IORB
END IF
END DO
C
C we may reorder the orbitals now
DO IORB=1,NBAS
DO IALPH=1,NBAS
CIN(IALPH,IORB)=CI(IALPH,ISRTI(IORB))
END DO
END DO
DO IORB=1,NBAS
DO IALPH=1,NBAS
CI(IALPH,IORB)=CIN(IALPH,IORB)
END DO
END DO
C
C we have to reorder as well IOCC and IOTYP
C this is done in the same array, we add 10*new_value and divide later
C by 10
C
C for NUMGRP we do the same with a factor of 1000
C
DO IORB=1,NBAS
IOCC(IORB)= IOCC(IORB)+10*MOD(IOCC(ISRTI(IORB)),10)
IOTYP(IORB)= IOTYP(IORB)+10*MOD(IOTYP(ISRTI(IORB)),10)
NUMGRP(IORB)= NUMGRP(IORB)+1000*MOD(NUMGRP(ISRTI(IORB)),1000)
DO J=1,4
IATYP(J,IORB)=IATYP(J,IORB)+1000*MOD(IATYP(J,ISRTI(IORB)),1000)
END DO
END DO
DO IORB=1,NBAS
IOCC(IORB)= IOCC(IORB)/10
IOTYP(IORB)= IOTYP(IORB)/10
NUMGRP(IORB)= NUMGRP(IORB)/1000
DO J=1,4
IATYP(J,IORB)=IATYP(J,IORB)/1000
END DO
END DO
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
FMO(IALPH,IBETA)=FAO(IALPH,IBETA)
END DO
END DO
CALL TRANS(FMO,CI,NBAS,NBAS)
DO I=1,NBAS
ORBEN(I)=FMO(I,I)
END DO
C rewrite the result of all the operations
WRITE(6,9400)
DO I=1,NBAS
IF (IOTYP(I).NE.5) THEN
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
$ ,(IATYP(J,I),J=1,IOTYP(I))
ELSE
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
END IF
END DO
C
C and now we treat the groups by the defined method, LCAN, LPIPM or LBOYS
C
IF (LCAN.OR.LPIPM.OR.LBOYS.OR.LEXTR) THEN
C//SKIPCP
cP IF (LEXTR) THEN
cP CALL EXTREM
cP CALL TIMING('EXTR')
cP ELSE
cP/ENDCP
DO IGRP=1,NGRP
NORBG=IFIGR(IGRP)-ISTGR(IGRP)+1
CD WRITE(6,*) ' NORBG = ',NORBG
IF (NORBG.NE.0) THEN
IF (LCAN) THEN
CALL GENCAN(FAO,FMO,CI(1,ISTGR(IGRP)),ORBEN(ISTGR(IGRP))
- ,NORBG)
ELSE
IF (LPIPM) THEN
CALL PIPEKM(CI(1,ISTGR(IGRP)),NORBG)
ELSE
IF (LBOYS) THEN
CALL BOYS(CI(1,ISTGR(IGRP)),NORBG)
ELSE
IF (LEXTR) THEN
CALL EXTRP(CI(1,ISTGR(IGRP)),NORBG)
END IF
END IF
END IF
END IF
END IF
END DO
C//SKIPCP
cP END IF
cP/ENDCP
ELSE
WRITE(6,*) ' no method specified, storing '
$ ,'orbitals in the obtained order'
END IF
C
C write the vector after ordering with respect to IOCC
CALL VSRTO(CI,IOCC,IOCCS)
C normalize the orbitals
CALL SNORM(CI,NBAS,NBAS)
C
C//SKIPCP
cP IF (.NOT.LEXTR) THEN
cP/ENDCP
DO I=1,NBAS
DO J=1,NBAS
FMO(I,J)=FAO(I,J)
END DO
END DO
CALL TRANS(FMO,CI,NBAS,NBAS)
DO I=1,NBAS
ORBEN(I)=FMO(I,I)
END DO
CALL BUBBLE(NBAS,NBAS,ORBEN,FMO,CI)
WRITE(6,*) ' new diagonal elements of F: '
WRITE(6,'(6(I5,F12.6))') (I,ORBEN(I),I=1,NBAS)
C
SOCC=0.D0
DO I=1,NOCC
SOCC=SOCC+ORBEN(I)
END DO
SVIRT=0.D0
DO I=1+NOCC,NBAS
SVIRT=SVIRT+ORBEN(I)
END DO
WRITE(6,*)
WRITE(6,*) ' sum of orbital energies: occupied ',SOCC
WRITE(6,*) ' virtuals ',SVIRT
WRITE(6,*)
C//SKIPCP
cP END IF
cP/ENDCP
C
C write vector
IF (LCAN ) CALL WVEC('CANL')
IF (LBOYS) CALL WVEC('BOYS')
IF (LPIPM) CALL WVEC('PIPM')
IF (LEXTR) CALL WVEC('EXTO')
C
ELSE
WRITE(6,*)
WRITE(6,*) ' no groups defined, leaving orbitals unchanged'
WRITE(6,*)
C
C get the dipole matrix elements
C
call baryzcal(CI,NBAS,baryz)
C
C we evaluate the Boys measure:
CALL BOYSME(NBAS,baryz)
C
END IF
C
IF (LSTAT) THEN
IF (LCAN ) THEN
CALL STATAN('CANL')
ELSE
IF (LBOYS) THEN
CALL STATAN('BOYS')
ELSE
IF (LPIPM) THEN
CALL STATAN('PIPM')
ELSE
IF (LEXTR) THEN
CALL STATAN('EXTO')
ELSE
CALL STATAN('XXXX')
END IF
END IF
END IF
END IF
END IF
IF (L6D) THEN
WRITE(6,*) ' projection of the occupied orbitals '
$ ,'onto cartesian Gaussians '
CALL PROJ6D
END IF
IF (LQMC.AND..NOT.LEXTR) CALL OUTQMC
CALL SPREAD(CI)
C rewrite the result of all the operations
WRITE(6,9400)
DO I=1,NBAS
IF (IOTYP(I).NE.5) THEN
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
$ ,(IATYP(J,I),J=1,IOTYP(I))
ELSE
WRITE(6,9401) I,ORBEN(I),IOCC(I),NUMGRP(I),COTYP(IOTYP(I))
END IF
END DO
WRITE(6,*)
WRITE(6,*)
C
IF (LMOLD) THEN
CALL FFCAL
IF (LBOYS) THEN
CALL MOLDEN('BOYS')
ELSE IF (LPIPM) THEN
CALL MOLDEN('PIPM')
ELSE IF (LEXTR) THEN
CALL MOLDEN('EXTO')
ELSE
CALL MOLDEN('CANL')
END IF
END IF
C
X=CPTIME(4)
CALL DATING(PNAME,2)
STOP
901 CONTINUE
WRITE(6,*) ' NO FILE FOUND'
STOP
902 CONTINUE
WRITE(6,*) ' NO FILE FOUND'
STOP
6100 CONTINUE
WRITE(6,*) ' NO FILE as starting vector FOUND ... '
STOP
END
C
SUBROUTINE ORTHO(CI,NOCC,NBAS)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(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),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION CI(NBASM,NVECT)
C
C normalisation of vectors
C
WRITE(6,*) ' norming the orbitals '
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,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(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,*)
C
DO I=1,NBAS
DO J=1,NBAS
CWRK(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(CWRK,CI,NOCC,NBAS)
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),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
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 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$$$ SUBROUTINE DENS(CI)
c$$$ INCLUDE 'param.h'
c$$$ INCLUDE 'common_syst.h'
c$$$ INCLUDE 'common_intu.h'
c$$$ INCLUDE 'flow.h'
c$$$ DIMENSION CI(NBASM,NBASM)
c$$$C
c$$$C calculate number of electrons
c$$$C
c$$$ WRITE(6,*) 'BASIS FUNCTION ASSIGNED NUMBER OF ELECTRONS'
c$$$ FTOT=0.D0
c$$$ DO K=1,NBAS
c$$$ FF=0.D0
c$$$ DO L=1,NBAS
c$$$ FF=FF+P(K,L)*SAO(K,L)
c$$$ END DO
c$$$ FF=FF
c$$$ FTOT=FTOT+FF
c$$$ IF (ABS(FF).GT.1.D-6) WRITE(6,'(5X,I6,15X,F13.6)') K,FF
c$$$ END DO
c$$$ WRITE(6,*)
c$$$ WRITE(6,*) 'TOTAL NUMBER OF ELECTRONS',FTOT
c$$$ WRITE(6,*)
c$$$ RETURN
c$$$ END
C
SUBROUTINE PFUNC(CI,NVECT,NBAS)
INCLUDE 'param.h'
DIMENSION CI(NBASM,NBASM)
C
DO IVEC=1,NVECT
WRITE(6,9901) IVEC
WRITE(6,9902) (IALPH,CI(IALPH,IVEC),IALPH=1,NBAS)
END DO
9901 FORMAT(' Function :',I5)
9902 FORMAT(5(I5,F10.5))
RETURN
END
C
SUBROUTINE PFUNCL(CI,NVECT,NBAS)
INCLUDE 'param.h'
DIMENSION CI(NBASM,NBASM)
C
DO IVEC=1,NVECT
WRITE(6,9901) IVEC
WRITE(6,9902) (IALPH,CI(IALPH,IVEC),IALPH=1,NBAS)
END DO
9901 FORMAT('# Function :',I5)
9902 FORMAT(I5,F23.17)
RETURN
END
C
SUBROUTINE RDMAT(NAME,ARRAY)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.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
IUNITR=31
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 RDINP
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(NBASM)
C.. INCLUDE 'common_mezey.h'
CHARACTER*4 KEYW(2),STR3
CHARACTER*6 KEYOPT(15),STR6
CHARACTER*80 LINE
DIMENSION INGRP(NBASM)
DATA KEYW /'*LPB','*END'/
DATA KEYOPT /'GROUP ','CANONI','PIPEKM','BOYS ','MOLDEN',
- 'FRAGME','EXTREM','QMC ','STATIS','NOVIRT',
- 'WIDTHG','CUTORB','6 D ','SMILES','XXXXXX'/
C
C DEFAULTS
LCAN=.FALSE.
LPIPM=.FALSE.
LBOYS=.FALSE.
LMOLD=.FALSE.
LFRAGM=.FALSE.
LEXTR=.FALSE.
LQMC=.FALSE.
NGRP=0
LSTAT=.FALSE.
LNOVIRT=.FALSE.
L6D=.FALSE.
CUTORB=1.D-8
LSMILE=.FALSE.
DO I=1,NBASM
NUMGRP(I)=0
END DO
WIDTHG=0.5D0
C
C structure of input:
C keyword for each program
IOINP=83
OPEN(UNIT=IOINP,FILE='INPUT.LPB',ERR=2217,FORM='FORMATTED',
- STATUS='OLD')
WRITE(6,*) ' found file INPUT.LPB'
C first, look for keyword 'LPB'
1 CONTINUE
READ(IOINP,'(A4)',END=2,ERR=921) STR3
IF (STR3.EQ.KEYW(1)) THEN
C look for Keyoptions
11 CONTINUE
READ(IOINP,'(A6)',END=2,ERR=921) STR6
WRITE(6,*) ' READ LINE |',STR6,'|'
IF (STR6(1:4).EQ.KEYW(2)) GO TO 2
IF (STR6.EQ.KEYOPT(1)) THEN
READ(IOINP,*) NGRP
WRITE(6,*) ' trying to read ',NGRP,' groups'
DO IGRP=1,NGRP
READ(IOINP,*) NGRP2
WRITE(6,*) ' trying to read ',NGRP2,' elements for group',IGRP
READ(IOINP,*) (INGRP(I),I=1,NGRP2)
DO I=1,NGRP2
NUMGRP(INGRP(I))=IGRP
END DO
END DO
ELSE IF (STR6.EQ.KEYOPT(2)) THEN
LCAN=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(3)) THEN
LPIPM=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(4)) THEN
LBOYS=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(5)) THEN
LMOLD=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(6)) THEN
LFRAGM=.TRUE.
READ(IOINP,*) NATDUM
READ(IOINP,*) (IFRAGM(I),I=1,NATDUM)
C extremely localized orbitals
ELSE IF (STR6.EQ.KEYOPT(7)) THEN
LEXTR=.TRUE.
C the QMC option
ELSE IF (STR6.EQ.KEYOPT(8)) THEN
LQMC=.TRUE.
C the statistical analysis
ELSE IF (STR6.EQ.KEYOPT(9)) THEN
LSTAT=.TRUE.
C no virtual orbitals
ELSE IF (STR6.EQ.KEYOPT(10)) THEN
LNOVIRT=.TRUE.
C width for Gaussian cut
ELSE IF (STR6.EQ.KEYOPT(11)) THEN
READ(IOINP,*) WIDTHG
C cutoff parameter for orbitals
ELSE IF (STR6.EQ.KEYOPT(12)) THEN
READ(IOINP,*) CUTORB
C projection on cartesian functions (1 S, 3 P, 6 D, 10 F, 15 G)
ELSE IF (STR6.EQ.KEYOPT(13)) THEN
L6D=.TRUE.
C reordering of oprbitals according to SMILES
ELSE IF (STR6.EQ.KEYOPT(14)) THEN
LSMILE=.TRUE.
ELSE
WRITE(6,*)
WRITE(6,*) ' UNKNOWN OPTION: |',STR6,'|'
WRITE(6,*)
WRITE(6,*) ' POSSIBLE OPTIONS IN THE BLOCK *LPB ... *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)
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'
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)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.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(FILEN)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
CHARACTER*4 FILEN
C
C output to a file
C
WRITE(6,*)
IUNITX=27
C
C localized orbitals by ORTHO
WRITE(6,*) ' WRITING VECTOR TO FILE '
OPEN(UNIT=IUNITX,FILE='VECTOR.'//FILEN,STATUS='UNKNOWN',FORM
$ ='FORMATTED')
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)
RETURN
END
C
SUBROUTINE WVCUT(FILEN)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
CHARACTER*4 FILEN
C
C output to a file
C
WRITE(6,*)
IUNITX=27
C
C localized orbitals by ORTHO
WRITE(6,*) ' WRITING VECTOR TO FILE '
OPEN(UNIT=IUNITX,FILE='VECTOR.'//FILEN//'.CUT',STATUS='UNKNOWN'
$ ,FORM='FORMATTED')
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)
RETURN
END
C
SUBROUTINE WVEXTR
C write extremely localized vector without orthogonalization
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
C
C output to a file
C
WRITE(6,*)
IUNITX=27
C
C localized orbitals by ORTHO
WRITE(6,*) ' WRITING VECTOR TO FILE '
OPEN(UNIT=IUNITX,FILE='VECTOR.EXTREME',STATUS='UNKNOWN'
-,FORM='FORMATTED')
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)
RETURN
END
C
SUBROUTINE GENCAN(FAO,FMO,CI,ORBEN,NVECT)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
DIMENSION FAO(NBASM,NBASM),FMO(NBASM,NBASM)
DIMENSION CI(NBASM,NVECT),ORBEN(NVECT)
DIMENSION CIN(NBASM,NVECT),CDUM1(NBASM),CDUM2(NBASM)
C
C here we diagonalize the Fock matrix, and store the new orbitals
C
DO I=1,NBAS
DO J=1,NBAS
FMO(I,J)=FAO(I,J)
END DO
END DO
CALL TRANS(FMO,CI,NVECT,NBAS)
C
IONE=1
CALL RS(NBASM,NVECT,FMO,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,NVECT)
WRITE(6,*)
C
C well, it should come down to a simple matrix multiplication ...
C
DO IALPH=1,NBAS
DO I=1,NVECT
CDUM=0.D0
DO J=1,NVECT
C the J components of eigenvector I
CDUM=CDUM+CIN(J,I)*CI(IALPH,J)
END DO
FMO(IALPH,I)=CDUM
END DO
END DO
C
DO IALPH=1,NBAS
DO I=1,NVECT
CI(IALPH,I)=FMO(IALPH,I)
END DO
END DO
C
RETURN
END
C
SUBROUTINE PIPEKM(CI,NVECT)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(NBASM)
C.. INCLUDE 'common_mezey.h'
DIMENSION CI(NBASM,NVECT)
LOGICAL LSWEEP
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
IF (NVECT.LE.1) THEN
WRITE(6,*) ' FOR NVECT=1 THERE IS NOTHING TO OPTIMIZE '
WRITE(6,*)
RETURN
END IF
CALL MEASUR(D,CI,NVECT)
WRITE(6,*) ' BEFORE: LOCALIZATION CRITERION D=',D
ITER=0
200 CONTINUE
ITER=ITER+1
CALL SWEEP(CI,NVECT,NBAS,NATOM,LSWEEP)
DOLD=D
CALL MEASUR(D,CI,NVECT)
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
RETURN
END
C
SUBROUTINE MEASUR(D,CI,NVECT)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(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.5D0*DBLE(NVECT)/D
RETURN
END
C
SUBROUTINE SWEEP(CI,NVECT,NBAS,NATOM,LSWEEP)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(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,NVECT
DO IT=IS+1,NVECT
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
S=DSIN(GAMMA)
S2=S*S
C2=1.D0-S2
C=DSQRT(C2)
SINA=S
COSA=C
c$$$ SINA=SIN(GAMMA)
c$$$ 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(NVECT*(NVECT-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.2.D-9) LSWEEP=.FALSE.
C
RETURN
END
C
SUBROUTINE SWEEPB(CI,NVECT,NBAS,LSWEEP)
INCLUDE 'param.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C.. INCLUDE 'common_boys.h'
DIMENSION CI(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,NVECT
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 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,NVECT
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(CI,NVECT)
C
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(NBASM)
C.. INCLUDE 'common_mezey.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C.. INCLUDE 'common_boys.h'
DIMENSION CI(NBASM,NVECT)
DIMENSION DIPDIA(NBASM,3)
LOGICAL LSWEEP
DIMENSION HDUM(NBASM,NBASM),ORD(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,*)
C
IF (NVECT.LE.1) THEN
WRITE(6,*) ' FOR NVECT=1 THERE IS NOTHING TO OPTIMIZE '
WRITE(6,*)
RETURN
END IF
IPRI=1
CALL MEASUR(D,CI,NVECT)
call baryzcal(CI,NVECT,DIPDIA)
C
C we evaluate the Boys measure:
CALL BOYSME(NVECT,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(CI,NVECT,NBAS,LSWEEP)
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 K=1,3
DO I=1,NVECT
DIPDIA(I,K)=DIPOL(I,I,K)
END DO
END DO
CALL BOYSME(NVECT,DIPDIA)
C
D=0.D0
CALL MEASUR(D,CI,NVECT)
C
RETURN
END
C
SUBROUTINE BOYSME(NORB,DIPDIA)
INCLUDE 'param.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C.. INCLUDE 'common_boys.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
C we have to have the diagonal of the dipole integrals in molecular orbitals
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
c$$$
c$$$ WRITE(6,9910)
c$$$ 9910 FORMAT(/,' barycenters of the orbitals (simple) '
c$$$ - ,//,5x,'No',10x,'x',15x,'y',15x,'z',/)
c$$$
c$$$ DO I=1,NORB
c$$$ xi=0.d0
c$$$ yi=0.d0
c$$$ zi=0.d0
c$$$ do ialph=1,nbas
c$$$ xalph=pos(1,iatz(ialph))
c$$$ yalph=pos(2,iatz(ialph))
c$$$ zalph=pos(3,iatz(ialph))
c$$$ c2=ci(ialph,i)*ci(ialph,i)
c$$$ xi=xi+c2*xalph
c$$$ yi=yi+c2*yalph
c$$$ zi=zi+c2*zalph
c$$$ end do
c$$$ WRITE(6,9902) I,xi,yi,zi
c$$$ END DO
c$$$
BBB=0.D0
DO K=1,3
DO I=1,NORB
DO J=I+1,NORB
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
C
SUBROUTINE SPREAD(CI)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
PARAMETER (NSECTM=30)
COMMON /MEZEY/ ISTRT(NATMX),NABAS(NATMX),NGRP,NUMGRP(NBASM)
$ ,ISTGR(NBASM),IFIGR(NBASM),IOTYP(NBASM),IATYP(4,NBASM)
$ ,IFRAGM(NATMX),ISMILE(NBASM)
C.. INCLUDE 'common_mezey.h'
DIMENSION CI(NBASM,NBAS),OCCP(NATMX),P(NBASM,NBASM),OCCT(NATMX)
C
DO IAT=1,NATOM
OCCT(IAT)=0.D0
END DO
WRITE(6,9900)(I,I=1,NATOM)
DO IORB=1,NBAS
C density matrix for each orbital
DO IATOM=1,NATOM
OCDUM=0.D0
DO IALPH=ISTRT(IATOM),ISTRT(IATOM)+NABAS(IATOM)
CIA=CI(IALPH,IORB)
DO IBETA=1,NBAS
OCDUM=OCDUM+CIA*CI(IBETA,IORB)*SAO(IALPH,IBETA)
END DO
END DO
OCCP(IATOM)=OCDUM
END DO
WRITE(6,9901) IORB,(2.D0*OCCP(I),I=1,NATOM)
IF (IORB.LE.NOCC) THEN
DO IAT=1,NATOM
OCCT(IAT)= OCCT(IAT) + 2.D0*OCCP(IAT)
END DO
END IF
C
C well, what type of orbital do we have?
OMX=0.D0
OMS=0.D0
DO IATOM=1,NATOM
OMS=OMS+OCCP(IATOM)
IF (OCCP(IATOM).GT.OMX) THEN
OMX=OCCP(IATOM)
IATYP(1,IORB)=IATOM
END IF
END DO
C how many other occupations are at least 50% of OMX?
IOMX=1
DO IATOM=1,NATOM
IF (OCCP(IATOM).GT.0.5D0*OMX.AND.IATOM.NE.IATYP(1,IORB)) THEN
IOMX=IOMX+1
IF (IOMX.LE.4) IATYP(IOMX,IORB)=IATOM
END IF
END DO
IOTYP(IORB)=MIN(5,IOMX)
IF (IOTYP(IORB).EQ.5) THEN
DO I=1,4
IATYP(I,IORB)=0
END DO
END IF
END DO
WRITE(6,*)
WRITE(6,9900) (IAT,IAT=1,NATOM)
WRITE(6,9902) (OCCT(IAT),IAT=1,NATOM)
9900 FORMAT(' Atom',I4,19(I8))
9901 FORMAT(I4,20(F8.4))
9902 FORMAT(' total ',20(F8.4))
C
RETURN
END
C
SUBROUTINE MOLDEN(FILEN)
C
C we try to create a MOLDEN file
C
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.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)
CHARACTER*4 FILEN
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
WRITE(6,*)
WRITE(6,*) ' writing file ',FILEN,'.MOLDEN '
WRITE(6,*)
C
IOUQMC=57
OPEN(UNIT=IOUQMC,FILE=FILEN//'.MOLDEN',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
CLOSE(IOUQMC,STATUS='DELETE')
OPEN(UNIT=IOUQMC,FILE=FILEN//'.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
ANTOAU=0.529177247D0
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
SUBROUTINE NORMBF(NSHLF,IPST)
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.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
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 EXTREM
C
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /EXTRE/ CIEX(NBASM,NBASM),CITMP(NBASM,NBASM)
C.. INCLUDE 'common_extreme.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION XLAMBD(NBASM,NBASM)
C
WRITE(6,*)
WRITE(6,*) ' creating projected occupied orbitals '
WRITE(6,*)
WRITE(6,*) ' the original orbitals '
CALL PFUNC(CI,NOCC,NBAS)
C XLAMBD is a working array, the first NOCC vectors are the new orbitals
OPEN(UNIT=25,FILE='PROJECTION.OCC',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
CALL PROJEC(CI,NOCC,NBAS,XLAMBD)
CLOSE(UNIT=25)
DO I=1,NOCC
ORBEN(I)=DBLE(I)
DO IALPH=1,NBAS
CI(IALPH,I)=XLAMBD(IALPH,I)
END DO
END DO
C invert notion occupied and virtual
IF (.NOT.LNOVIRT) THEN
WRITE(6,*) ' creating projected virtual orbitals '
WRITE(6,*)
NVIRT=NBAS-NOCC
OPEN(UNIT=25,FILE='PROJECTION.VIRT',FORM='FORMATTED',STATUS
$ ='UNKNOWN')
CALL PROJEC(CI(1,NOCC+1),NVIRT,NBAS,XLAMBD)
CLOSE(UNIT=25)
C
DO I=NOCC+1,NBAS
ORBEN(I)=DBLE(I)
DO IALPH=1,NBAS
CI(IALPH,I)=XLAMBD(IALPH,I-NOCC)
END DO
END DO
END IF
C
CALL PFUNC(CI,NBAS,NBAS)
C
CALL WVEXTR
IF (LSTAT) CALL STATAN('EXTN')
C
IF (LMOLD) THEN
CALL FFCAL
CALL MOLDEN('EXTN')
END IF
C
IF (LQMC) THEN
CALL OUTQMC
ELSE
C we orthogonalize and verify if the vector spaces were well orthogonal
WRITE(6,*) ' orthogonalizing occupied space '
CALL SHALF(CI,NOCC,NBAS)
C
NVIRT=NBAS-NOCC
WRITE(6,*) ' orthogonalizing virtual space '
CALL SHALF(CI(1,NOCC+1),NVIRT,NBAS)
C we use again the dummy array XLAMBD
DO IALPH=1,NBAS
DO IBETA=1,NBAS
XLAMBD(IALPH,IBETA)=SAO(IALPH,IBETA)
END DO
END DO
CALL TRANS(XLAMBD,CI,NBAS,NBAS)
C
C largest element coupling the two spaces
SMX=0.D0
DO I=1,NOCC
DO J=NOCC+1,NBAS
SSS=ABS(XLAMBD(I,J))
IF (SSS.GT.SMX) SMX=SSS
END DO
END DO
WRITE(6,*)
WRITE(6,*) ' largest overlap coupling occupied '
$ ,'and virtual orbital spaces: ',SMX
WRITE(6,*)
C
C and what about the Fock matrix elements?
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
FMO(IALPH,IBETA)=FAO(IALPH,IBETA)
END DO
END DO
C
CALL TRANS(FMO,CI,NBAS,NBAS)
C
FMX=0.D0
DO I=1,NOCC
DO J=NOCC+1,NBAS
FFF=ABS(FMO(I,J))
IF (FFF.GT.FMX) FMX=FFF
END DO
END DO
WRITE(6,*)
WRITE(6,*) ' largest Fock matrix element coupling occupied '
$ ,'and virtual orbital spaces: ',FMX
WRITE(6,*)
C
C write Fock matrix
c$$$ OPEN(UNIT=17,FILE='MATRIX',FORM='FORMATTED',STATUS='UNKNOWN')
c$$$ DO I=1,NBAS
c$$$ DO J=1,NBAS
c$$$ WRITE(17,*) I,J,FMO(I,J)
c$$$ END DO
c$$$ END DO
c$$$ CLOSE(17)
C
END IF
C
RETURN
END
C
SUBROUTINE PROJEC(CI,NOCC,NBAS,XLAMBD)
INCLUDE 'param.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /EXTRE/ CIEX(NBASM,NBASM),CITMP(NBASM,NBASM)
C.. INCLUDE 'common_extreme.h'
C
C on input we have CI(NBAS,*), on output XLAMBD, without having
C destroyed CI
C
DIMENSION XLAMBD(NBASM,NBASM),CI(NBASM,*)
C
C local arrays
C
DIMENSION XSPRD(NBASM),IPERM(NBASM),PTMP(NBASM,NBASM),CDUM1(NBASM)
$ ,CDUM2(NBASM)
CHARACTER*3 CPF
CPF=' | '
C
C we have to construct
C
C \Lambda_{\alpha\beta} = \delta_|alpha\beta}
C -\sum_{\gamma}S_{\alpha\gamma} P_{\beta\gamma}
C
C with P_{\beta\gamma} = \sum_{n\in virt} c_{\beta\,n} c_{\gamma\,n}
C
C in fact we may just think one second: P(virt) = 1 - P(occ)
C so we do not have to look at the virtuals for projecting onto the occupied
C
DO IALPH=1,NBAS
DO IBETA=1,NBAS
PTMP(IALPH,IBETA)=0.D0
XLAMBD(IALPH,IBETA)=0.D0
END DO
END DO
C
DO N=1,NOCC
DO IALPH=1,NBAS
CAN=CI(IALPH,N)
DO IBETA=1,NBAS
PTMP(IALPH,IBETA)=PTMP(IALPH,IBETA)+CAN*CI(IBETA,N)
END DO
END DO
END DO
C
DO IALPH=1,NBAS
DO IGAMM=1,NBAS
SSS=SAO(IALPH,IGAMM)
DO IBETA=1,NBAS
XLAMBD(IBETA,IALPH)= XLAMBD(IBETA,IALPH)+SSS*PTMP(IBETA,IGAMM)
END DO
END DO
END DO
C
write(6,*)
write(6,*) ' Projection: '
C
c$$$ CALL PFUNC(XLAMBD,NBAS,NBAS)
C
C XLAMBD(IALPH,IBETA) is not XLAMBD(IBETA,IALPH)
C
C the new orbitals are XLAMBD(IALPH,I)
C
C we may eliminate the projections with zero length, assuming that the
C AOs are not linearly dependent
C
NPROJ=NBAS
IVEC=1
100 CONTINUE
C calculate norm of vector IVEC
SSS=0.D0
DO IALPH=1,NBAS
SSS=SSS+XLAMBD(IALPH,IVEC)*XLAMBD(IALPH,IVEC)
END DO
IF (SSS.LT.1.D-12) THEN
C exchange with vector NPROJ
WRITE(6,*) ' VECTOR No',IVEC,' has zero length '
DO IALPH=1,NBAS
XLAMBD(IALPH,IVEC)=XLAMBD(IALPH,NPROJ)
END DO
NPROJ=NPROJ-1
C test again
IF (IVEC.LT.NPROJ) GO TO 100
END IF
IVEC=IVEC+1
IF (IVEC.LE.NPROJ) GO TO 100
C
C we end up here if IVEC=NPROJ
C
c$$$ CALL PFUNC(XLAMBD,NPROJ,NBAS)
WRITE(6,*)
WRITE(6,*) ' we have ',NPROJ,' remaining non-zero projections '
WRITE(6,*)
C
C here we have to implement the localization criterion
C
C a criterion: we divide by the largest, and sum the coefficients
C this will be largest for 1+1 (delocalized)
C and smallest for 1+0 (completely localized)
C
DO I=1,NPROJ
C find the largest
XMX=0.D0
DO IALPH=1,NBAS
IF (ABS(XLAMBD(IALPH,I)).GT.XMX) XMX=ABS(XLAMBD(IALPH,I))
END DO
C
XSP=0.D0
DO IALPH=1,NBAS
XLAMBD(IALPH,I)=XLAMBD(IALPH,I)/XMX
XSP=XSP+ABS(XLAMBD(IALPH,I))
c$$$ XSP=XSP+ABS(XLAMBD(IALPH,I))**0.25D0
END DO
XSPRD(I)=XSP
END DO
C
C we write the projections to file 25, already open
WRITE(25,*) NPROJ,NBAS
DO I=1,NPROJ
WRITE(25,9543) I
9543 FORMAT('# Function No ',I5)
WRITE(25,'(E21.12)') (XLAMBD(IALPH,I),IALPH=1,NBAS)
WRITE(25,*)
END DO
c$$$ CALL PFUNCL(XLAMBD,NPROJ,NBAS)
C
write(6,*) ' normalizing the projections '
CALL SNORM(XLAMBD,NPROJ,NBAS)
c$$$C here we have to implement the localization criterion
c$$$ DO I=1,NPROJ
c$$$ XSP=0.D0
c$$$ DO IALPH=1,NBAS
c$$$ XSP=XSP+XLAMBD(IALPH,I)**4.D0
c$$$ END DO
c$$$ XSPRD(I)=XSP
c$$$ END DO
C
C select the NOCC most localized ones
C order the orbitals according to the spread (heapsort)
C
CALL HPSORT(XSPRD,IPERM,NPROJ)
C and store them in the right order in CITMP
WRITE(6,*) ' Spread of the projected orbitals '
WRITE(6,*)
WRITE(6,'(4(I4,E12.4))') (I,XSPRD(I),I=1,NPROJ)
WRITE(6,*) ' Permutations: '
WRITE(6,*)
DO I=1,NPROJ,18
WRITE(6,'(1H|,18(I3,1H|))') (INDX,INDX=I,MIN(I+17,NPROJ))
WRITE(6,'(1H|,18(A3,1H|))') (CPF,INDX=I,MIN(I+17,NPROJ))
WRITE(6,'(1H|,18(I3,1H|))') (IPERM(INDX),INDX=I,MIN(I+17,NPROJ))
WRITE(6,*)
END DO
C
DO I=1,NPROJ
DO IALPH=1,NBAS
c$$$ CITMP(IALPH,IPERM(I))=XLAMBD(IALPH,I)
CITMP(IALPH,I)=XLAMBD(IALPH,IPERM(I))
END DO
END DO
C
C recalculate spread ....
C
DO I=1,NPROJ
C find the largest
XMX=0.D0
DO IALPH=1,NBAS
IF (ABS(CITMP(IALPH,I)).GT.XMX) XMX=ABS(CITMP(IALPH,I))
END DO
XSP=0.D0
C number of very small coefficients
NSM=0
DO IALPH=1,NBAS
XSP=XSP+ABS(CITMP(IALPH,I))/XMX
c$$$ XSP=XSP+ABS(XLAMBD(IALPH,I))**0.25D0
IF (ABS(CITMP(IALPH,I)).LT.1.D-6) NSM=NSM+1
END DO
write(6,*) ' spread of function ',I,': ',XSP,NSM
END DO
C
CALL PFUNC(CITMP,NPROJ,NBAS)
C
C in PTMP is the overlap matrix of the 'new' orbitals
C
C now we start with the first two vectors and add vectors as long as we
C do not have NOCC different vectors. After adding we check if the set
C is still not linearly dependent. If yes, we skip the last vector added
C and take the next one.
C
C a first threshold, if this is too tight, we may lower it and restart
C the whole procedure
C
THRESH=1.D-6
2000 CONTINUE
NVEC=2
IOFF=0
1002 CONTINUE
C Calculate the determinant of the new overlap matrix
DO I=1,NBAS
DO J=1,NBAS
PTMP(I,J)=SAO(I,J)
END DO
END DO
CALL TRANS(PTMP,CITMP,NVEC,NBAS)
C
C we need only the eigenvalues, and only the smallest of the symmetric
C matrix
C the sorting array XSRT is not necessary any more, we use it for the eigenvalues
C XLAMBD becomes a dummy array
IMATZ=0
CALL RS(NBASM,NVEC,PTMP,XSPRD,IMATZ,XLAMBD,CDUM1,CDUM2,IERR)
C
WRITE(6,*) ' eigenvalues: '
WRITE(6,'(4(I4,E12.4))') (I,XSPRD(I),I=1,NVEC)
WRITE(6,*)
IF (XSPRD(1).LE.THRESH) THEN
IOFF=IOFF+1
IF (NVEC+IOFF.GT.NBAS) THEN
C if we ran out of vectors we have to loosen the threshold and restart
THRESH=THRESH*0.1D0
WRITE(6,*)
WRITE(6,*) ' restarting with a new threshold: ',thresh
WRITE(6,*)
GO TO 2000
END IF
WRITE(6,*) ' adding function No ',NVEC+IOFF,' for ',NVEC
WRITE(6,*)
DO IALPH=1,NBAS
CDUM=CITMP(IALPH,NVEC)
CITMP(IALPH,NVEC)=CITMP(IALPH,NVEC+IOFF)
CITMP(IALPH,NVEC+IOFF)=CDUM
END DO
GO TO 1002
ELSE
NVEC=NVEC+1
IF (NVEC.LE.NOCC) THEN
WRITE(6,*) ' adding function No ',NVEC
GO TO 1002
END IF
END IF
C
c$$$ WRITE(6,*)
c$$$ WRITE(6,*)
c$$$ WRITE(6,*) ' after the reduction: '
c$$$ WRITE(6,*)
c$$$ WRITE(6,*)
c$$$ CALL PFUNC(CITMP,NOCC,NBAS)
C
C normalize them for conserving the Slater determinant
C the result is returned in XLAMBD
C
C assigning functions, we use again XLAMBD for that purpose
C
C we don''t have to change the functions, it is just sufficient to
C assure that there is a 1-to-1 assignement for the normation
C
DO I=1,NOCC
DO J=1,NOCC
XLAMBD(I,J)=0.D0
END DO
END DO
C
DO I=1,NOCC
DO IALPH=1,NBAS
CAI=CI(IALPH,I)
DO IBETA=1,NBAS
SAI=CAI*SAO(IALPH,IBETA)
DO J=1,NOCC
CBJ=CITMP(IBETA,J)
XLAMBD(I,J)=XLAMBD(I,J)+SAI*CBJ
END DO
END DO
END DO
END DO
C
WRITE(6,*)
WRITE(6,*) ' the overlap matrix: horiz.: '
$ ,'read vectors, vert.: projections '
WRITE(6,*)
CALL PFUNC(XLAMBD,NOCC,NOCC)
WRITE(6,*)
C
DO I=1,NOCC
INDX=0
INDXO=0
SSS=0.D0
DO J=I,NOCC
IF (ABS(XLAMBD(J,I)).GT.SSS) THEN
INDX=J
SSS=ABS(XLAMBD(J,I))
END IF
END DO
9215 FORMAT(' exchanging row ',I4,' and ',I5,': elements ',2F14.9)
C is the element negative?
IF (XLAMBD(INDX,I).LE.0.D0) THEN
MULT=-1
ELSE
MULT=1
END IF
XMULT=DBLE(MULT)
IF (INDX.NE.I) THEN
C exchange row I and row INDX in HERMAT
WRITE(6,9215) I,MULT*INDX,XLAMBD(I,I),XLAMBD(INDX,I)
DO K=1,NOCC
CDUM=XLAMBD(I,K)
XLAMBD(I,K)=XLAMBD(INDX,K)*XMULT
XLAMBD(INDX,K)=CDUM
END DO
ELSE
IF (MULT.EQ.-1) THEN
D=-D
WRITE(6,9215) I,-I,XLAMBD(I,I),-XLAMBD(I,I)
DO K=1,NOCC
XLAMBD(I,K)=XMULT*XLAMBD(I,K)
END DO
END IF
END IF
END DO
C
WRITE(6,*)
WRITE(6,*)
$ ' diagonal of the overlap matrix after reordering '
WRITE(6,'(5(I4,F14.8))') (I,XLAMBD(I,I),I=1,NOCC)
WRITE(6,*)
C
DO I=1,NOCC
SSS=XLAMBD(I,I)
DO IALPH=1,NBAS
XLAMBD(IALPH,I)=CITMP(IALPH,I)/SSS
END DO
END DO
C
CALL TIMING('PROJ')
C
RETURN
END
C
SUBROUTINE HPSORT(XSRT,IWRK,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C sorting after criterion XSRT corresponding array of dimension N
C IWRK contains the permutation upon completion
DIMENSION XSRT(N),IWRK(N)
C
DO I=1,N
IWRK(I)=I
END DO
C
C SORT BY XSRT
IF (N.LT.2) RETURN
L=N/2+1
IR=N
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
IRA=IWRK(L)
RRA=XSRT(L)
ELSE
IRA=IWRK(IR)
RRA=XSRT(IR)
XSRT(IR)=XSRT(1)
IWRK(IR)=IWRK(1)
IR=IR-1
IF (IR.EQ.1) THEN
IWRK(1)=IRA
XSRT(1)=RRA
RETURN
END IF
END IF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF(J.LT.IR)THEN
IF(XSRT(J).LT.XSRT(J+1))J=J+1
END IF
IF(RRA.LT.XSRT(J))THEN
IWRK(I)=IWRK(J)
XSRT(I)=XSRT(J)
I=J
J=J+J
ELSE
J=IR+1
END IF
GOTO 20
END IF
XSRT(I)=RRA
IWRK(I)=IRA
GOTO 10
C
END
C
SUBROUTINE OUTQMC
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.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(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)
DO IKK=-ILL+1,ILL-1
IBAS=IBAS+1
WRITE(6,9221) IBAS,IAT,CST(ILL),IKK
END DO
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
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
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
END DO
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) EHF,ECORR,EHF+ECORR,PHP
8009 FORMAT(' Hartree-Fock Energy: ',/,F20.12,/
$ ,' Correlation Energy by the selected determinants: ',/,F20
$ .12,/,' Total Energy: ',/,F20.12,/,' /: '
$ ,/,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
WRITE(IINQMC,8111) INUM,INUM+NJASIM,IZERO,ONE
9000 FORMAT(' Determinant ',I5,' Weight and coefficients ',/,F20.12)
C
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
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)
WRITE(IINQMC,8112) IORB,'alpha spin',(OWBUF(IPRIM)
$ ,IPRIM=1,IPCOUNT)
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)
WRITE(IINQMC,8112) IORB,'beta spin',(OWBUF(IPRIM)
$ ,IPRIM=1,IPCOUNT)
END DO
9901 FORMAT(' Molecular orbital No ',I3,' ',A,/,3(F20.12))
8112 FORMAT(' Molecular orbital No ',I3,' ',A,/,3(F23.17))
CLOSE(IOUQMC)
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 baryzcal(CI,NVECT,BARYZ)
INCLUDE 'param.h'
COMMON /CBOYS/ DIPOL(NBASM,NBASM,3)
C.. INCLUDE 'common_boys.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
DIMENSION CI(NBASM,NVECT),BARYZ(NBASM,3)
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,NVECT,NBAS)
CALL TRANS(DIPOL(1,1,2),CI,NVECT,NBAS)
CALL TRANS(DIPOL(1,1,3),CI,NVECT,NBAS)
DO K=1,3
DO I=1,NVECT
BARYZ(I,K)=DIPOL(I,I,K)
END DO
END DO
WRITE(6,*)
WRITE(6,*)
WRITE(6,9909)
9909 FORMAT(/,' baryzcal: barycenters of the orbitals '
- ,//,5x,'No',10x,'x',15x,'y',15x,'z',/)
DO I=1,NVECT
WRITE(6,9902) I,(baryz(I,K),K=1,3)
9902 FORMAT(I6,3F16.8)
END DO
WRITE(6,*)
WRITE(6,*)
return
end
SUBROUTINE STATAN(type)
INCLUDE 'param.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
LOGICAL LCUTCI(NBASM,NBASM)
character*4 type
C we have to analyze the occupied and the virtual orbitals separately
C save the orbitals on file
open(unit=12,file='vtmp',status='unknown',form='unformatted')
write(12) ((ci(ialph,i),ialph=1,nbas),i=1,nbas)
close(12)
C divide by the largest
do i=1,nbas
xmx=0.d0
do ialph=1,nbas
if (abs(ci(ialph,i)).gt.xmx) xmx=abs(ci(ialph,i))
end do
do ialph=1,nbas
ci(ialph,i)=ci(ialph,i)/xmx
end do
end do
WRITE(6,*)
WRITE(6,*) ' simple statistical analysis '
WRITE(6,*) ' - - - - - - - - - - - - - - - - - - - '
WRITE(6,*)
WRITE(6,*)
C
call statdetail
C
write(6,*)
write(6,*) ' statistical analysis including distances '
write(6,*)
C renormalize the orbitals
call snorm(ci,nbas,nbas)
C calculate the barycenters of the orbitals
CALL BARYZCAL(ci,nbas,baryz)
CUTlarge=1.d-2
C a model exponent
expcut=widthg*widthg
write(6,*)
write(6,*) ' length unit is ',1.d0/widthg,' a.u.'
write(6,*) ' Gaussian exp(-r^2/(2 sigma)) with sigma = ',0.5d0
$ /expcut
write(6,*)
C we scale every coefficient by exp(-expcut*(r-r_i)^2)
C but only if it is smaller than 10-2
do i=1,nbas
bx=baryz(i,1)
by=baryz(i,2)
bz=baryz(i,3)
C write(6,*) ' orbital i: bx by bz: ',i,bx,by,bz
do ialph=1,nbas
LCUTCI(IALPH,I)=.FALSE.
if (abs(ci(ialph,i)).lt.cutlarge) then
rx=pos(1,iatz(ialph))
ry=pos(2,iatz(ialph))
rz=pos(3,iatz(ialph))
rr= (rx-bx)*(rx-bx) + (ry-by)*(ry-by) + (rz-bz)*(rz-bz)
rex=expcut*rr
if (rex.gt.25.d0) then
ci(ialph,i)=0.d0
else
ci(ialph,i)=ci(ialph,i)*exp(-expcut*rr)
end if
end if
end do
end do
do i=1,nbas
do ialph=1,nbas
if (abs(ci(ialph,i)).lt.cutorb) LCUTCI(IALPH,I)=.TRUE.
end do
end do
C divide by the largest
do i=1,nbas
xmx=0.d0
do ialph=1,nbas
if (abs(ci(ialph,i)).gt.xmx) xmx=abs(ci(ialph,i))
end do
do ialph=1,nbas
ci(ialph,i)=ci(ialph,i)/xmx
end do
end do
call statdetail
C restore the orbitals
open(unit=12,file='vtmp',status='old',form='unformatted')
read(12) ((ci(ialph,i),ialph=1,nbas),i=1,nbas)
close(12)
C put zeroes for the tails cut
do i=1,nbas
do ialph=1,nbas
if (lcutci(ialph,i)) ci(ialph,i)=0.d0
end do
end do
CALL WVCUT(TYPE)
C restore the orbitals
open(unit=12,file='vtmp',status='old',form='unformatted')
read(12) ((ci(ialph,i),ialph=1,nbas),i=1,nbas)
close(12,status='delete')
end
subroutine statdetail
INCLUDE 'param.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
PARAMETER (NCASE=20)
DIMENSION KOCC(NCASE,NBASM)
do i=1,nbas
do icase=1,ncase
kocc(icase,i)=0
end do
do ialph=1,nbas
call getcas(ci(ialph,i),icase,ncase)
c$$$ write(6,*) ialph,i,icase,orbital(ialph,i)
kocc(icase,i)=kocc(icase,i)+1
end do
end do
izero=0
write(6,'(i3,i7,19i4)') izero,(icase,icase=1,ncase)
write(6,*) ' occupied orbitals '
write(6,'(i4,i7,19i4)') (i,(kocc(icase,i),icase=1,ncase),i=1,nocc)
do i=2,nocc
do icase=1,ncase
kocc(icase,1)=kocc(icase,1)+kocc(icase,i)
end do
end do
write(6,*) ' sums '
C write(6,'(i4,i7,19i5)') izero,(kocc(icase,1),icase=1,ncase)
write(6,'(i7,i7)') (icase, kocc(icase,1),icase=1,ncase)
write(6,*) ' '
do icase=2,ncase
kocc(icase,1)=kocc(icase-1,1)+kocc(icase,1)
end do
write(6,*) ' cumulatif '
C write(6,'(i4,i7,19i6)') izero,(kocc(icase,1),icase=1,ncase)
write(6,'(2i7)') (icase,kocc(icase,1),icase=1,ncase)
write(6,*) ' '
write(6,*) ' '
IF (.not.lnovirt) then
write(6,*) ' virtual orbitals '
write(6,'(i4,i7,19i4)') (i,(kocc(icase,i),icase=1,ncase),i=nocc+1
$ ,nbas)
do icase=1,ncase
kocc(icase,1)=0
end do
do i=1+nocc,nbas
do icase=1,ncase
kocc(icase,1)=kocc(icase,1)+kocc(icase,i)
end do
end do
write(6,*) ' sums '
write(6,'(i4,i7,19i6)') izero,(kocc(icase,1),icase=1,ncase)
write(6,*) ' '
do icase=2,ncase
kocc(icase,1)=kocc(icase-1,1)+kocc(icase,1)
end do
write(6,*) ' cumulatif '
write(6,'(i4,i7,19i7)') izero,(kocc(icase,1),icase=1,ncase)
write(6,*) ' '
write(6,*) ' '
end if
return
end
subroutine getcas(x,i,nmax)
implicit double precision (a-h,o-z)
amin=10**dble(-nmax)
if (abs(x).lt.amin) then
i=nmax
else
if (abs(x).gt.0.1d0) then
i=1
else
xl10=log(10.D0)
xl=log(abs(x))/xl10
il=1-xl
i=min(il,nmax)
i=max(i,1)
end if
end if
return
end
C
SUBROUTINE EXTRP(CIV,NVECT)
C
INCLUDE 'param.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /EXTRE/ CIEX(NBASM,NBASM),CITMP(NBASM,NBASM)
C.. INCLUDE 'common_extreme.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /INTU/ SAO(NBASM,NBASM),FAO(NBASM,NBASM),FMO(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
DIMENSION XLAMBD(NBASM,NBASM)
DIMENSION CIV(NBASM,NVECT)
C
WRITE(6,*)
WRITE(6,*) ' projected orbitals '
WRITE(6,*)
C XLAMBD is a working array, the first NOCC vectors are the new orbitals
CALL PROJEC(CIV,NVECT,NBAS,XLAMBD)
C the projection is in XLAMBD
DO IVECT=1,NVECT
DO IALPH=1,NBAS
CI(IALPH,IVECT)=XLAMBD(IALPH,IVECT)
END DO
END DO
C
RETURN
END
C
SUBROUTINE PROJ6D
INCLUDE 'param.h'
C
COMMON /FLOW/ EPS,WIDTHG,CUTORB,LCAN,LPIPM,LMOLD,LBOYS,LFRAGM
$ ,LEXTR,LQMC,LSTAT,LNOVIRT,L6D,LSMILE
LOGICAL LCAN,LPIPM,LMOLD,LBOYS,LFRAGM,LEXTR,LQMC,LSTAT
$ ,LNOVIRT,L6D,LSMILE
C.. INCLUDE 'flow.h'
COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT,POS(3,NATMX),NZ(NATMX)
$ ,IATZ(NBASM)
COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),NSH(NATMX)
$ ,NPRIM(NSHLMX),IL(NSHLMX),NPX(NLMAX,NATMX)
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),BARYZ(NBASM,3),IOCC(NBASM)
$ ,IOCCS(NBASM)
C.. INCLUDE 'common_vect.h'
PARAMETER (NDIM6D=NBASM+20)
DIMENSION VECT6D(NDIM6D,NBASM),DENS6D(NDIM6D,NDIM6D)
C
C this should become the interface routine to the Maeder-Claverie
C program MTPCOR
WRITE(6,*)
WRITE(6,*) ' projection onto cartesian functions '
WRITE(6,*)
$ ' preparation of input to multipole'
$ ,' decomposition of Claverie'
WRITE(6,*)
C
IUNITW=13
OPEN(UNIT=IUNITW,FILE='INPUT.MTPCOR',STATUS='UNKNOWN'
$ ,FORM='FORMATTED')
WRITE(IUNITW,9100) 'INTINP'
9100 FORMAT(' &',A6,/,' NEWSTYLE=T ',/,' &END')
C
C we add the SYSTEM.ORTHO FILE
IPRIM=0
NBAS6D=0
NSHL=0
WRITE(IUNITW,*) NATOM
DO IAT=1,NATOM
WRITE(IUNITW,9101) NZ(IAT),NSH(IAT),(POS(J,IAT),J=1,3)
DO ISHL=1,NSH(IAT)
NSHL=NSHL+1
ILL=IL(NSHL)
IF (ILL.EQ.0) THEN
NBAS6D=NBAS6D+1
ELSE IF (ILL.EQ.1) THEN
NBAS6D=NBAS6D+3
ELSE IF (ILL.EQ.2) THEN
NBAS6D=NBAS6D+6
ELSE IF (ILL.EQ.3) THEN
NBAS6D=NBAS6D+10
ELSE IF (ILL.EQ.4) THEN
NBAS6D=NBAS6D+15
ELSE
WRITE(6,*) ' not yet beyond g functions '
END IF
WRITE(IUNITW,'(2I6)') ILL,NPRIM(NSHL)
DO III=0,NPRIM(NSHL)-1
IPRIM=IPRIM+1
WRITE(IUNITW,'(F24.12,F20.12)') EXX(IPRIM),COEFF(IPRIM)
END DO
END DO
END DO
9101 FORMAT(I3,I6,3F20.12)
C
C write a dummy number for NORB in mtpcor
C
WRITE(13,*) 123456
C
C write second namelist for mtpcor
C
WRITE(13,9100) 'PRPINP'
C
WRITE(6,*) ' we have ',NBAS ,' spherical basis functions '
WRITE(6,*) ' we have ',NBAS6D,' cartesian basis functions '
C
C do the expansion from spherical to cartesian Gaussians
C
NSHL=0
NBB=0
NBB6D=0
DO IAT=1,NATOM
DO ISHL=1,NSH(IAT)
NSHL=NSHL+1
ILL=IL(NSHL)
C S
IF (ILL.EQ.0) THEN
NBAS6D=NBAS6D+1
NBB=NBB+1
NBB6D=NBB6D+1
DO IORB=1,NOCC
VECT6D(NBB6D,IORB)=CI(NBB,IORB)
END DO
C P
ELSE IF (ILL.EQ.1) THEN
C px, py, pz
DO III=1,3
NBB=NBB+1
NBB6D=NBB6D+1
DO IORB=1,NOCC
VECT6D(NBB6D,IORB)=CI(NBB,IORB)
END DO
END DO
C D
ELSE IF (ILL.EQ.2) THEN
C xy, yz, 3z2-r2, xz, x2-y2
C
C xx, yy, zz, xy, xz, yz
C
DO IORB=1,NOCC
VECT6D(NBB6D+1,IORB)=CI(NBB+5,IORB)-CI(NBB+3,IORB)
VECT6D(NBB6D+2,IORB)=-CI(NBB+5,IORB)-CI(NBB+3,IORB)
VECT6D(NBB6D+3,IORB)=2.D0*CI(NBB+1,IORB)
VECT6D(NBB6D+4,IORB)=CI(NBB+1,IORB)
VECT6D(NBB6D+5,IORB)=CI(NBB+4,IORB)
VECT6D(NBB6D+6,IORB)=CI(NBB+2,IORB)
END DO
NBB=NBB+5
NBB6D=NBB6D+6
ELSE
WRITE(6,*) ' projection not yet beyond f functions '
END IF
END DO
END DO
c$$$
c$$$ DALTON output
c$$$
c$$$ +----------------------------------------------+
c$$$ ! Cartesian to spherical transformation matrix !
c$$$ +----------------------------------------------+
c$$$
c$$$ Moment order: 2
c$$$
c$$$ Column 1 Column 2 Column 3 Column 4
c$$$ 1 0.00000000 0.00000000 -1.00000000 0.00000000
c$$$ 2 1.00000000 0.00000000 0.00000000 0.00000000
c$$$ 3 0.00000000 0.00000000 0.00000000 1.00000000
c$$$ 4 0.00000000 0.00000000 -1.00000000 0.00000000
c$$$ 5 0.00000000 1.00000000 0.00000000 0.00000000
c$$$ 6 0.00000000 0.00000000 2.00000000 0.00000000
c$$$
c$$$ Column 5
c$$$ 1 1.00000000
c$$$ 4 -1.00000000
c$$$
c$$$
c$$$ +----------------------------------------------+
c$$$ ! Cartesian to spherical transformation matrix !
c$$$ +----------------------------------------------+
c$$$
c$$$ Moment order: 3
c$$$
c$$$ Column 1 Column 2 Column 3 Column 4
c$$$ 2 3.00000000 0.00000000 -1.00000000 0.00000000
c$$$ 3 0.00000000 0.00000000 0.00000000 1.00000000
c$$$ 5 0.00000000 1.00000000 0.00000000 0.00000000
c$$$ 7 -1.00000000 0.00000000 -1.00000000 0.00000000
c$$$ 8 0.00000000 0.00000000 0.00000000 1.00000000
c$$$ 9 0.00000000 0.00000000 4.00000000 0.00000000
c$$$ 10 0.00000000 0.00000000 0.00000000 -0.66666667
c$$$
c$$$
c$$$
c$$$ +----------------------------------------------+
c$$$ ! Cartesian to spherical transformation matrix !
c$$$ +----------------------------------------------+
c$$$
c$$$ Moment order: 4
c$$$
c$$$ Column 1 Column 2 Column 3 Column 4
c$$$ 2 1.000000E+00 0.000000E+00 -1.000000E+00 0.000000E+00
c$$$ 5 0.000000E+00 3.000000E+00 0.000000E+00 -1.000000E+00
c$$$ 7 -1.000000E+00 0.000000E+00 -1.000000E+00 0.000000E+00
c$$$ 9 0.000000E+00 0.000000E+00 6.000000E+00 0.000000E+00
c$$$ 12 0.000000E+00 -1.000000E+00 0.000000E+00 -1.000000E+00
c$$$ 14 0.000000E+00 0.000000E+00 0.000000E+00 1.333333E+00
c$$$
c$$$ Column 5 Column 6 Column 7 Column 8
c$$$ 1 -1.250000E-01 0.000000E+00 -2.058292E+16 0.000000E+00
c$$$ 3 0.000000E+00 -1.000000E+00 0.000000E+00 -3.333333E-01
c$$$ 4 -2.500000E-01 0.000000E+00 -1.000000E+00 0.000000E+00
c$$$ 6 1.000000E+00 0.000000E+00 1.234975E+17 0.000000E+00
c$$$ 8 0.000000E+00 -1.000000E+00 0.000000E+00 1.000000E+00
c$$$ 10 0.000000E+00 1.333333E+00 0.000000E+00 0.000000E+00
c$$$ 11 -1.250000E-01 0.000000E+00 2.058292E+16 0.000000E+00
c$$$ 13 1.000000E+00 0.000000E+00 -1.234975E+17 0.000000E+00
c$$$ 15 -3.333333E-01 0.000000E+00 0.000000E+00 0.000000E+00
c$$$
c$$$ Column 9
c$$$ 1 -1.666667E-01
c$$$ 4 1.000000E+00
c$$$ 11 -1.666667E-01
C
C create density matrix
DO IALPH=1,NBB6D
DO IBETA=IALPH,NBB6D
DDD=0.D0
DO IORB=1,NOCC
DDD=DDD+VECT6D(IALPH,IORB)*VECT6d(IBETA,IORB)
END DO
DENS6D(IALPH,IBETA)=2.D0*DDD
END DO
END DO
WRITE(IUNITW,9102) ((DENS6D(IALPH,IBETA),IBETA=IALPH,NBB6D)
$ ,IALPH=1,NBB6D)
9102 FORMAT(4E20.12)
C
CLOSE(IUNITW)
C
RETURN
END