#define _MCSCF_ 0 #define _RHF_ 1 #define _UHF_ 2 #define _ROHF_ 3 #define _ROHFSTD_ 4 #define _ROHFSCN_ 5 *********************************************************************** module scf_embed *********************************************************************** * Contains (almost) all SCF-related embedding procedures * 2023-05-12 Bence Hégely *********************************************************************** implicit none contains subroutine readh_embedding_alphabeta $(Hcore_embed,Hcore_embed_beta,r8heap,ifile,symtra,nbasis,offset, $lsa) implicit none ! intent(in) integer, intent(in) :: ifile integer, intent(in) :: nbasis integer, intent(in) :: offset double precision, intent(in) :: r8heap(*) double precision, intent(in) :: symtra(nbasis,nbasis) logical, intent(in) :: lsa ! intent(out) double precision, intent(out) :: Hcore_embed(*) double precision, intent(out) :: Hcore_embed_beta(*) ! open(ifile,file='OEINT',form='UNFORMATTED',action='READ', $ position='REWIND',status='OLD') read(ifile) call readone $( Hcore_embed,ifile,symtra,nbasis,r8heap,offset,lsa) call readone $( Hcore_embed_beta,ifile,symtra,nbasis,r8heap,offset,lsa) close(ifile) end subroutine readh_embedding_alphabeta double precision function calc_embedding_uhf_ener $(dft,pcm,nuc,sqrsize,p_alpha,p_beta,fock_alpha,fock_beta, $hcore_alpha,hcore_beta,exc,ddot) result( ener ) implicit none character(len=*) :: dft character(len=*) :: pcm integer, intent(in) :: sqrsize double precision, intent(in) :: nuc double precision, intent(in) :: exc double precision, intent(in) :: p_alpha(*) double precision, intent(in) :: p_beta(*) double precision, intent(in) :: fock_alpha(*) double precision, intent(in) :: fock_beta(*) double precision, intent(in) :: hcore_alpha(*) double precision, intent(in) :: hcore_beta(*) double precision, external :: ddot if(trim(dft).eq.'off'.and.trim(pcm).eq.'off') then ener = nuc $ + 0.5d0*ddot( sqrsize , p_alpha , 1 , fock_alpha , 1 ) $ + 0.5d0*ddot( sqrsize , p_beta , 1 , fock_beta , 1 ) $ + ddot( sqrsize , p_alpha , 1 , hcore_alpha , 1 ) $ + ddot( sqrsize , p_beta , 1 , hcore_beta , 1 ) else ener = nuc $ + exc $ + ddot( sqrsize , p_alpha , 1 , hcore_alpha , 1 ) $ + ddot( sqrsize , p_beta , 1 , hcore_beta , 1 ) endif end function calc_embedding_uhf_ener C subroutine projection_embed_ROUTE_em2 $(SCFTYPE,ifltln,iintln,iout,scrfile1,oeintfile,varsfile,offset, $nbasis,nal,nbe,nelec,nfroz,sqrsize,embed,dual,core,pcm,dft,route, $moselect,subnuc,ener,exc,symtra,p_subsys,p_subsys_beta, $fock_subsys,fock_subsys_beta,r8heap,r8heap1n,r8heap2n,r8heap3n, $r8heap_rs,sorig,lsa,ldual,lembed_grad,itol) c use common_mod, only: ifltln,iintln,iout,scrfile1,oeintfile, c & varsfile implicit none integer, intent(in) :: SCFTYPE integer, intent(in) :: ifltln integer, intent(in) :: iintln integer, intent(in) :: iout integer, intent(in) :: scrfile1 integer, intent(in) :: oeintfile integer, intent(in) :: varsfile integer, intent(in) :: offset integer, intent(in) :: nbasis integer, intent(in) :: nal integer, intent(in) :: nbe integer, intent(in) :: nelec integer, intent(in) :: nfroz integer, intent(in) :: sqrsize character(len=*), intent(in) :: embed character(len=*), intent(in) :: dual character(len=*), intent(in) :: core character(len=*), intent(in) :: pcm character(len=*), intent(in) :: dft character(len=*), intent(in) :: route character(len=*), intent(in) :: moselect double precision, intent(in) :: subnuc double precision, intent(in) :: ener double precision, intent(in) :: exc double precision, intent(in) :: itol double precision :: symtra(nbasis,nbasis) double precision, intent(in) :: sorig(*) double precision :: p_subsys(*) double precision :: p_subsys_beta(*) double precision :: fock_subsys(*) double precision :: fock_subsys_beta(*) double precision, intent(in), target :: r8heap(sqrsize) double precision, intent(in), target :: r8heap1n(sqrsize) double precision, intent(in), target :: r8heap2n(sqrsize) double precision, intent(in), target :: r8heap3n(sqrsize) double precision, intent(in), target :: r8heap_rs(sqrsize) logical, intent(in) :: lsa logical, intent(in) :: ldual logical, intent(in) :: lembed_grad integer :: ncore double precision :: ener_subsys_low double precision :: ENER_subsysAlow_and_corr double precision, external :: ddot double precision :: ENER_supersys_low double precision, pointer :: Hcore_subsys(:)=> null() double precision, pointer :: fock_embed(:) => null() double precision, pointer :: fock_embed_beta(:) $ => null() ! call get_subsys_oeint call calc_subsys_energy_lowlevel call save_f2_a ! only for gradient calculation Hcore_subsys => null() call print_subsys_energy_lowlevel call save_subsys_energy_lowlevel call build_and_save_embedding_potential call overwrite_oeint call update_varsfile call calc_first_order_correction call save_first_order_correction ! contains subroutine get_subsys_oeint Hcore_subsys => r8heap open(oeintfile,file='OEINT',form='UNFORMATTED', $ position='REWIND',status='OLD',action='READ') read(oeintfile) if(moselect.eq.'ecore') then ! Actually, this is the supersystem h_core call roeint $ (r8heap3n,r8heap3n,Hcore_subsys,oeintfile,nbasis) else read(oeintfile) call roeint $ (r8heap3n,r8heap3n,Hcore_subsys,oeintfile,nbasis) endif close(oeintfile) end subroutine get_subsys_oeint ! subroutine calc_subsys_energy_lowlevel if( SCFTYPE .eq. _RHF_ ) then if(trim(dft).eq.'off'.and.trim(pcm).eq.'off') then ener_subsys_low $ = subnuc $ + 0.5d0 * ddot $ (sqrsize , p_subsys , 1 , fock_subsys ,1 ) $ + ddot $ (sqrsize , p_subsys , 1 , Hcore_subsys ,1 ) else ener_subsys_low $ = subnuc $ + exc $ + ddot $ (sqrsize , p_subsys , 1 , Hcore_subsys ,1 ) endif else if( SCFTYPE .eq. _UHF_ ) then if(trim(dft).eq.'off'.and.trim(pcm).eq.'off') then ener_subsys_low $ = subnuc $ + 0.5d0*ddot $ (sqrsize , p_subsys , 1 , fock_subsys ,1 ) $ + 0.5d0*ddot $ (sqrsize , p_subsys_beta , 1 , fock_subsys_beta,1 ) $ + ddot $ (sqrsize , p_subsys , 1 , Hcore_subsys ,1 ) $ + ddot $ (sqrsize , p_subsys_beta , 1 , Hcore_subsys ,1 ) else ener_subsys_low $ = subnuc $ + exc $ + ddot $ (sqrsize , p_subsys , 1 , Hcore_subsys ,1 ) $ + ddot $ (sqrsize , p_subsys_beta , 1 , Hcore_subsys ,1 ) endif else write(iout,'(a)') $' Unknown SCF type for embedding calculation.' call mrccend(1) endif end subroutine calc_subsys_energy_lowlevel ! subroutine print_subsys_energy_lowlevel write(iout,*) write(iout, &"(' ***ENERGY OF THE ACTIVE SUBSYSTEM ',F28.16,' [AU]')") & ener_subsys_low write(iout,'(a)') $' (low level) ' end subroutine print_subsys_energy_lowlevel ! subroutine save_subsys_energy_lowlevel open(varsfile,file='VARS',form='UNFORMATTED', $ position='APPEND',status='OLD',action='WRITE') write(varsfile) 'edftalow ',ifltln,ener_subsys_low close(varsfile) end subroutine save_subsys_energy_lowlevel ! subroutine build_and_save_embedding_potential ! Read low level AO Fock of the supersystem double precision, pointer :: fock_supersys_low(:) $ =>null() double precision, pointer :: fock_supersys_low_beta(:) $ =>null() fock_supersys_low => r8heap open(scrfile1,file='FOCK',form='UNFORMATTED', $ status='OLD',position='REWIND',action='READ') call roeint $ (r8heap3n,r8heap3n,fock_supersys_low,scrfile1,nbasis) if( SCFTYPE .eq. _UHF_ ) then fock_supersys_low_beta => r8heap1n call roeint $(r8heap3n,r8heap3n,fock_supersys_low_beta,scrfile1,nbasis) endif if( ldual ) then call getvar('ncorenew ',ncore) else read(scrfile1) ncore if(embed.eq.'huzinaga') then if(core.eq.'corr ') ncore=0 ncore=ncore+nfroz ! warning, bug detector beeps endif endif close(scrfile1) fock_embed => r8heap call daxpy(sqrsize,-1.d0,fock_subsys,1,fock_embed,1) if( SCFTYPE .eq. _UHF_ ) then fock_embed_beta => r8heap1n call daxpy $ (sqrsize,-1.d0,fock_subsys_beta,1,fock_embed_beta,1) fock_supersys_low_beta => null() endif fock_supersys_low => null() end subroutine build_and_save_embedding_potential ! subroutine overwrite_oeint open(oeintfile,file='OEINT_EMBED',form='UNFORMATTED', $ status='REPLACE',action='WRITE',position='REWIND') call woeint $ (r8heap3n,r8heap3n,sorig,oeintfile,0.0d0,nbasis) call woeint $ (r8heap3n,r8heap3n,fock_embed,oeintfile,0.0d0,nbasis) if( SCFTYPE .eq. _UHF_ ) then call woeint $ (r8heap3n,r8heap3n,fock_embed_beta,oeintfile,0.0d0,nbasis) endif close(oeintfile) end subroutine ! subroutine update_varsfile if(embed.ne.'fdm ') then call ishell('cp VARS VARS.old') open(varsfile,file='VARS.new',form='unformatted') if(embed.ne.'huzinaga') then write(varsfile) 'nal ',iintln,nal write(varsfile) 'nbe ',iintln,nbe write(varsfile) 'nelec ',iintln,nelec endif write(varsfile) 'ncore ',iintln,ncore close(varsfile) call ishell('cat VARS.new VARS.old > VARS') call ishell('rm -f VARS.new VARS.old') endif end subroutine update_varsfile ! subroutine calc_first_order_correction double precision, pointer :: Hcore_supersys(:) =>null() double precision, pointer :: p_subsys_improved(:) =>null() Hcore_supersys => r8heap2n call readh $ (Hcore_supersys,oeintfile,symtra,nbasis,offset,lsa) call daxpy(sqrsize,-1.d0,Hcore_supersys,1,fock_embed,1) if( SCFTYPE .eq. _UHF_ ) call daxpy $ (sqrsize,-1.d0,Hcore_supersys,1,fock_embed_beta,1) call getvar('edftab ', ENER_supersys_low) if( SCFTYPE .eq. _RHF_ ) then ENER_subsysAlow_and_corr $ = ENER_supersys_low $ - ddot(sqrsize,p_subsys,1,fock_embed,1) $ - ener else if( SCFTYPE .eq. _UHF_ ) then ENER_subsysAlow_and_corr $ = ENER_supersys_low $ - ddot(sqrsize,p_subsys,1,fock_embed,1) $ - ddot(sqrsize,p_subsys_beta,1,fock_embed_beta,1) $ - ener endif if( ldual ) then p_subsys_improved => r8heap_rs call daxpy $ (sqrsize,1.d0,Hcore_supersys,1,fock_subsys,1) open(scrfile1,file='SCFDENSITIES',status='OLD', $ position='REWIND',form='UNFORMATTED') read(scrfile1) call roeint $ (r8heap3n,r8heap3n,p_subsys_improved,scrfile1,nbasis) close(scrfile1) call daxpy $ (sqrsize,-1.d0,p_subsys,1,p_subsys_improved,1) ENER_subsysAlow_and_corr = ENER_subsysAlow_and_corr $ - ddot(sqrsize,p_subsys_improved,1,fock_subsys,1) $ - ddot(sqrsize,p_subsys_improved,1,fock_embed,1) p_subsys_improved => null() endif call ishell('mv OEINT_EMBED OEINT') if( SCFTYPE .eq. _UHF_ ) then fock_embed_beta => null() endif fock_embed => null() Hcore_supersys => null() end subroutine calc_first_order_correction ! subroutine save_first_order_correction open(varsfile,file='VARS',form='UNFORMATTED', $ position='APPEND',status='OLD',action='WRITE') write(varsfile) $ 'edftbcor ',ifltln,ENER_subsysAlow_and_corr close(varsfile) end subroutine save_first_order_correction ! subroutine save_f2_a if(.not.lembed_grad) return * Save F_2[D^A] for Huzinaga embedding gradient open(scrfile1,file='FOCK2_A',form='unformatted') call dcopy(nbasis*nbasis,fock_subsys,1,r8heap1n,1) call daxpy $ (nbasis*nbasis,1.0d0,Hcore_subsys,1,r8heap1n,1) call woeint $ (r8heap2n,r8heap2n,r8heap1n,scrfile1,itol,nbasis) call woeint $ (r8heap2n,r8heap2n,Hcore_subsys,scrfile1,itol,nbasis) close(scrfile1) end subroutine save_f2_a end subroutine projection_embed_ROUTE_em2 ! subroutine print_energies_of_partitioned_orbitals $(iout,qmreg,lhuzi,embed,huzitype,eigenvalue,mapping, $nfroz,nvfroz,nvfrozb,nal,nbe,nbasis,ndao,ndocc,scftype,route) implicit none integer, intent(in) :: scftype integer, intent(in) :: iout integer, intent(in) :: huzitype integer, intent(in) :: nfroz integer, intent(in) :: nvfroz integer, intent(in) :: nvfrozb integer, intent(in) :: nal integer, intent(in) :: nbe integer, intent(in) :: nbasis integer, intent(in) :: ndao integer, intent(in) :: ndocc integer, intent(in) :: mapping(*) double precision, intent(in) :: eigenvalue(*) character(len=4), intent(in) :: route character(len=8), intent(in) :: embed character(len=8), intent(in) :: qmreg logical, intent(in) :: lhuzi character(len=32) :: label character(len=32) :: label_beta integer :: i,q q = merge( 2 , 1 , scftype .eq. _UHF_ ) write(iout,*) if( scftype .eq. _RHF_ ) then if(qmreg.eq.'0 '.or.lhuzi.or.embed.eq.'scl ') $ label = ' (frozen orbital) ' if(route.eq.'sch ') $ label = ' (orbital of the environment)' write(iout,*) 'Orbital energies [au]:' write(iout,*) ' Irrep Alpha Beta' do i=1,nbasis-ndao label = ' ' if( i .ge. 1 .and. i .le. nfroz ) then if(route.eq.'sch ') $ label = ' (orbital of the environment)' if(qmreg.eq.'0 '.or. $ lhuzi.or.embed.eq.'scl ') $ label = ' (frozen orbital) ' endif if( huzitype .eq. 2 $ .and. $ i .ge. nal + 1 .and. i .le. nal + nvfroz ) then label = ' (frozen virtual) ' endif write(*,"(i5,3x,a4,2f15.6,a32)") i,'A ',!irlab(mosym(mapping(I))), $EIGENVALUE(MAPPING(I)), $EIGENVALUE(MAPPING(I+(Q-1)*nbasis)+(Q-1)*nbasis),label if(i.eq.ndocc) write(iout,*) enddo else if ( scftype .eq. _UHF_ ) then write(iout,*) 'Orbital energies [au]:' write(iout,*) $' Irrep Alpha Beta' do i=1,nbasis-ndao label = ' ' label_beta = ' ' if( i .ge. 1 .and. i .le. nfroz .and. lhuzi ) then label = ' (frozen) ' label_beta = ' (frozen) ' endif if( huzitype .eq. 2 ) then if( i .gt. nal .and. i .le. nal + nvfroz ) then label = ' (frozen) ' endif if( i .gt. nbe .and. i .le. nbe + nvfrozb ) then label_beta = ' (frozen) ' endif endif write(*,"(i5,3x,a4,f15.6,a12,f15.6,a12)") i,'A ',!irlab(mosym(mapping(I))), $EIGENVALUE(MAPPING(I)),label, $EIGENVALUE(MAPPING(I+(Q-1)*nbasis)+(Q-1)*nbasis),label_beta if(i.eq.ndocc) write(iout,*) enddo endif end subroutine print_energies_of_partitioned_orbitals C subroutine build_huzinaga_operator $(scftype,huzitype,route,nbasis,sqrsize,fock,rs,r8heap) ************************************************************************ * Build Huzinaga operator * Note: if virtual orbitals are projected, than a slightly different * operator is built to match the requirements of the sign of * the usual virtuals. ************************************************************************ implicit none integer, intent(in) :: scftype integer, intent(in) :: huzitype character(len=*), intent(in) :: route integer, intent(in) :: nbasis integer, intent(in) :: sqrsize double precision, intent(inout) :: r8heap(*) double precision, intent(inout) :: fock(*) double precision, intent(in) :: rs(*) if( scftype .eq. _RHF_ ) then call huzproj(nbasis,fock,rs,r8heap) if( huzitype .eq. 2 ) then call huzproj(nbasis,fock,rs(sqrsize+1),r8heap) call huzproj(nbasis,fock,rs(sqrsize+1),r8heap) endif else if( scftype .eq. _UHF_ ) then call huzproj $( nbasis , fock , rs , r8heap ) call huzproj $( nbasis , fock(sqrsize+1) , rs(sqrsize+1) , r8heap ) if(huzitype.eq.2) then call huzproj $( nbasis , fock , rs(2*sqrsize+1) , r8heap ) call huzproj $( nbasis , fock , rs(2*sqrsize+1) , r8heap ) call huzproj $( nbasis , fock(sqrsize+1) , rs(3*sqrsize+1) , r8heap ) call huzproj $( nbasis , fock(sqrsize+1) , rs(3*sqrsize+1) , r8heap ) endif endif end subroutine build_huzinaga_operator ! subroutine huzproj(nbasis,fock,rs,work) implicit none integer, intent(in) ::nbasis double precision, intent(inout) :: fock(*) double precision, intent(in) :: rs(*) double precision, intent(in) :: work(*) call dcopy(nbasis**2,fock,1,work,1) call dsyr2k('u','t',nbasis,nbasis,-1.d0,rs,nbasis,work,nbasis, $1.d0,fock,nbasis) call filllo(fock,nbasis) end subroutine huzproj C ************************************************************************ * Read the [ projector ] X [AO overlap] for the Huzinaga equations ************************************************************************ subroutine read_huzproj_file $(r8heap,rs,scrfile1,nbf,scftype,sqrsize,huzitype) implicit none integer, intent(in) :: huzitype integer, intent(in) :: scrfile1 integer, intent(in) :: nbf(5) integer, intent(in) :: scftype integer, intent(in) :: sqrsize double precision, intent(inout) :: rs(*) double precision, intent(in) :: r8heap(*) open(scrfile1,file='HUZPROJ',form='UNFORMATTED',action='READ', $ status='OLD',position='REWIND') call rtdmx(r8heap,r8heap,rs,scrfile1,nbf(5),nbf(5)) if( SCFTYPE .eq. _UHF_ ) $ call rtdmx $(r8heap,r8heap,rs(sqrsize+1),scrfile1,nbf(5),nbf(5)) if( huzitype .eq. 2 .and. SCFTYPE .eq. _RHF_ ) $ call rtdmx $(r8heap,r8heap,rs(sqrsize+1),scrfile1,nbf(5),nbf(5)) if( huzitype .eq. 2 .and. SCFTYPE .eq. _UHF_ ) $ call rtdmx $(r8heap,r8heap,rs(2*sqrsize+1),scrfile1,nbf(5),nbf(5)) if( huzitype .eq. 2 .and. SCFTYPE .eq. _UHF_ ) $ call rtdmx $(r8heap,r8heap,rs(3*sqrsize+1),scrfile1,nbf(5),nbf(5)) close(scrfile1) end subroutine read_huzproj_file C ************************************************************************ * Sort orbitals if the Huzinaga equation is used. * Note: the projected orbitals change their signs thus the * standard ordering from pseig mixes the frozen occupied and * virtual orbitals ************************************************************************ subroutine sort_orbitals_for_embedding(scftype, huzitype, $nbasis, sqrsize, c, rs, eigenvalue, r8heap, r8heap1n, iout, $nfroz, nvfroz, nvfrozb, sorig, qmreg, route, nal, nbe , $qmreg_froz) implicit none integer, intent(in) :: scftype integer, intent(in) :: huzitype integer, intent(in) :: nbasis integer, intent(in) :: sqrsize integer, intent(in) :: iout integer, intent(in) :: nfroz integer, intent(in) :: nvfroz integer, intent(in) :: nvfrozb integer, intent(in) :: nal, nbe integer, intent(in) :: qmreg_froz character(len=4), intent(in) :: route character(len=8), intent(in) :: qmreg double precision, intent(in) :: rs(*) double precision, intent(in) :: sorig(*) double precision, intent(in) :: r8heap(*) double precision, intent(in) :: r8heap1n(*) double precision, intent(inout) :: eigenvalue(*) double precision, intent(inout) :: c(*) ! sort occupied orbitals (rhf case) / sort alpha occupied orbitals (uhf case) call huzord $(nbasis, c, rs, eigenvalue, r8heap, r8heap1n, iout, nfroz, sorig, $qmreg, route, nal, qmreg_froz, .false. ) if( scftype .eq. _UHF_ ) then ! sort beta occupied orbitals call huzord $(nbasis, c(sqrsize+1), rs(sqrsize+1), eigenvalue(nbasis+1), $r8heap, r8heap1n,iout, nfroz, sorig, qmreg, route, nbe,0, .false.) endif if( huzitype.eq.2 ) then if ( scftype .eq. _RHF_ ) then ! sort alpha virtual orbitals (rhf case) call huzord $( nbasis, c, rs(sqrsize+1), eigenvalue, r8heap, r8heap1n, iout, $nvfroz, sorig, qmreg, route, nal, 0, .true. ) else if ( scftype .eq. _UHF_ ) then ! sort alpha virtual orbitals (uhf case) call huzord $( nbasis, c, rs(2*sqrsize+1), eigenvalue, r8heap, r8heap1n, iout, $nvfroz, sorig, qmreg, route, nal, 0, .true.) ! sort beta virtual orbitals (uhf case) call huzord $( nbasis, c(sqrsize+1), rs(3*sqrsize+1), eigenvalue(nbasis+1), $r8heap, r8heap1n, iout, nvfrozb, sorig, qmreg, route, nbe, 0 $,.true.) endif endif end subroutine sort_orbitals_for_embedding ************************************************************************ subroutine huzord(nbasis,c,rs,eig,r8heap2a_nbfXnbf, $r8heap2b_nbfXnbf,iout,nfroz,s,qmreg,route,nelec,qmreg_froz, $sort_virtuals) ************************************************************************ * Move frozen orbitals to the beginning in the case of Huzinaga eqs ************************************************************************ implicit none ! intent(in) character(len=4), intent(in) :: route character(len=8), intent(in) :: qmreg integer, intent(in) :: nelec integer, intent(in) :: nbasis integer, intent(in) :: nfroz integer, intent(in) :: iout integer, intent(in) :: qmreg_froz double precision, intent(in) :: s(nbasis,nbasis) double precision, intent(in) :: rs(nbasis,nbasis) double precision, intent(in), target :: r8heap2a_nbfXnbf $ (nbasis,nbasis) double precision, intent(in), target :: r8heap2b_nbfXnbf $ (nbasis,nbasis) logical, intent(in) :: sort_virtuals ! intent(inout) double precision, intent(inout) :: eig(nbasis) double precision, intent(inout) :: c(nbasis,nbasis) ! scr double precision, pointer :: work(:,:) =>null() double precision, pointer :: c_envXs(:,:) =>null() double precision, pointer :: MO_overlap(:) =>null() double precision, pointer :: rscr_nbf(:) =>null() ! external double precision, external :: ddot ! local integer :: fst_ind, lst_ind,i call init_pointers call calc_MO_overlap if( .not. sort_virtuals ) then call sort_occupied_MOs else call sort_virtual_MOs endif call check_for_intruder_states call print_warnings call set_nullpointers C contains subroutine init_pointers work => r8heap2a_nbfXnbf c_envXs => r8heap2b_nbfXnbf MO_overlap => r8heap2a_nbfXnbf(:,1) rscr_nbf => r8heap2a_nbfXnbf(:,2) work = 0.0d0 end subroutine init_pointers subroutine set_nullpointers work => null() c_envXs => null() MO_overlap => null() rscr_nbf => null() end subroutine set_nullpointers subroutine calc_MO_overlap integer :: iMO call dgemm('n','n',nbasis,nbasis,nbasis,1.d0,rs,nbasis,c, $nbasis,0.d0,work,nbasis) call dsymm('l','u',nbasis,nbasis,1.d0,s,nbasis,work, $nbasis,0.d0,c_envXs,nbasis) MO_overlap = 0.0d0 rscr_nbf = 0.0d0 do iMO = 1 , nbasis MO_overlap( iMO ) = ddot $ ( nbasis ,c_envXs(1,iMO) ,1 ,c(1,iMO),1) enddo end subroutine calc_MO_overlap subroutine sort_occupied_MOs fst_ind = 1 lst_ind = nfroz call sort_orbitals( get_largest_MO_overlap_ind,(/1,nbasis/)) call sort_orbitals( get_highest_eig_ind, (/ 1 , lst_ind /) ) fst_ind = nfroz + 1 lst_ind = nbasis call sort_orbitals( get_lowest_eig_ind, (/ 1 , lst_ind /) ) end subroutine sort_occupied_MOs subroutine sort_virtual_MOs fst_ind = nelec + 1 lst_ind = nelec + nfroz call sort_orbitals(get_largest_MO_overlap_ind,(/1,nbasis /)) call sort_orbitals( get_lowest_eig_ind, (/ 1 , lst_ind /) ) fst_ind = nelec + nfroz + 1 lst_ind = nbasis call sort_orbitals( get_lowest_eig_ind, (/ 1 , lst_ind /) ) end subroutine sort_virtual_MOs subroutine sort_orbitals( func , func_range ) integer, external :: func integer, intent(in) :: func_range(2) integer :: i,ii double precision :: tmp i = fst_ind do while( i .le. lst_ind ) ii = func( i , func_range(2) ) if( i .eq. ii ) continue tmp=eig(i) eig(i)=eig(ii) eig(ii)=tmp rscr_nbf(1:nbasis)=c(1:nbasis,i) c(1:nbasis,i)=c(1:nbasis,ii) c(1:nbasis,ii)=rscr_nbf(1:nbasis) tmp=MO_overlap(i) MO_overlap(i)=MO_overlap(ii) MO_overlap(ii)=tmp i = i + 1 enddo end subroutine sort_orbitals pure integer function get_largest_MO_overlap_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = maxloc(MO_overlap(from:to),1)+from-1 end function get_largest_MO_overlap_ind pure integer function get_lowest_MO_overlap_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = minloc(MO_overlap(from:to),1)+from-1 end function get_lowest_MO_overlap_ind pure integer function get_highest_eig_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = maxloc(eig(from:to),1)+from-1 end function get_highest_eig_ind pure integer function get_lowest_eig_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = minloc(eig(from:to),1)+from-1 end function get_lowest_eig_ind subroutine check_for_intruder_states integer :: nMOs integer :: maxlap_ind, minlap_ind double precision :: tmp, maxlap, minlap integer :: i fst_ind = merge( nelec + 1 , 1 , sort_virtuals ) lst_ind = merge( nelec + nfroz , nfroz , sort_virtuals ) minlap = $ MO_overlap( get_lowest_MO_overlap_ind( fst_ind, lst_ind ) ) fst_ind = merge( nelec + nfroz + 1 ,nfroz + 1, sort_virtuals ) lst_ind = nbasis maxlap_ind = get_largest_MO_overlap_ind( fst_ind, lst_ind ) maxlap = MO_overlap( maxlap_ind ) nMOs = merge( nbasis - nelec , nelec , sort_virtuals ) if( nfroz .ne. 0 .and. nfroz .ne. nMOs ) then if( maxlap .gt. 0.1d0 ) then write(iout,'(a)') $' Warning! Intruder states are detected!' write(iout,'(a)') $' This can cause convergence problems.' write(iout,'(a,f10.5,a,i8,a)') $' The largest overlap amongst intruder states: ', $maxlap,' (MO: ',maxlap_ind,')' endif tmp = maxlap - minlap if( abs(tmp) .lt. 0.1d0 ) then write(iout,*) write(iout,'(a)') $' Warning! Non-projector orbitals critically'// $' overlap with the projector.' write(iout,'(a)') $' This can cause convergence problems.' write(iout,'(a,f10.5)') $' Difference of orbital overlaps = ',abs(tmp) write(iout,*) endif endif end subroutine check_for_intruder_states subroutine print_warnings integer :: pMO fst_ind = merge( nelec + 1 , 1 , sort_virtuals ) lst_ind = merge( nelec + nfroz + 1, nfroz , sort_virtuals ) if( qmreg .eq. '0' .and. .not.sort_virtuals $ .and. route .eq.'em3' ) then lst_ind = lst_ind - qmreg_froz endif do pMO = fst_ind , lst_ind if( eig( pMO ) .lt. 0.0d0 ) then call warning_msg( pMO ) endif enddo end subroutine print_warnings subroutine warning_msg( pMO ) integer, intent(in) :: pMO character(len=8) :: cscr write(cscr,'(i8)') pMO if( .not. sort_virtuals ) then write(iout,'(a,a,a)') $' Warning! The orbital energy of frozen occupied orbital ' $// trim(adjustl(cscr)) // ' is negative!' else write(iout,'(a)') $' Warning! The orbital energy of frozen virtual orbital ' $// trim(adjustl(cscr)) // ' is negative!' endif end subroutine warning_msg end subroutine huzord C subroutine build_operators_for_sch $(r8heap,fock,rs,p,sorig,nbasis,sqrsize) integer, intent(in) :: nbasis integer, intent(in) :: sqrsize double precision, intent(in) :: r8heap(*) double precision, intent(inout) :: fock(*) double precision, intent(in) :: p(*) double precision, intent(inout) :: rs(*) double precision, intent(in) :: sorig(*) ! fock : embedded subsystem A Fock matrix ! fock(sqrsize+1) : supersystem low-level Fock matrix ! p : supersystem density (AB) ! p(sqrsize+1) : subsystem density B ! rs : (projector of subsystem B)*S ! rs(sqrsize+1) : (projector of subsystem A)*S call dsymm('l','u',nbasis,nbasis,0.5d0,p(sqrsize+1),nbasis, $sorig,nbasis,0.d0,rs,nbasis) call huzproj(nbasis,fock,rs,r8heap) call dcopy(nbasis**2,p,1,r8heap,1) call daxpy(sqrsize,-1.d0,p(sqrsize+1),1,r8heap,1) call dsymm('l','u',nbasis,nbasis,0.5d0,r8heap,nbasis,sorig, $nbasis,0.d0,rs(sqrsize+1),nbasis) call huzproj(nbasis,fock(sqrsize+1),rs(sqrsize+1),r8heap) end subroutine build_operators_for_sch C subroutine write_RMS_of_commutators_for_sch(iout,devcomm) implicit none integer, intent(in) :: iout double precision, intent(in) :: devcomm(2) write(iout, &"(' RMS of [subsystem Fock, supersystem P] ',f25.14)") $devcomm(1) write(iout, &"(' RMS of [supersystem Fock, supersystem P] ',f25.14)") $devcomm(2) end subroutine write_RMS_of_commutators_for_sch C subroutine sort_orbitals_for_sch_embedding $(scftype, huzitype, $nbasis, sqrsize, c, rs, eigenvalue, r8heap, r8heap1n, iout, $nfroz, nvfroz, nvfrozb, sorig, qmreg, route, nal, nbe , $qmreg_froz) implicit none integer, intent(in) :: scftype integer, intent(in) :: huzitype integer, intent(in) :: nbasis integer, intent(in) :: sqrsize integer, intent(in) :: iout integer, intent(in) :: nfroz integer, intent(in) :: nvfroz integer, intent(in) :: nvfrozb integer, intent(in) :: nal, nbe integer, intent(in) :: qmreg_froz character(len=4), intent(in) :: route character(len=8), intent(in) :: qmreg double precision, intent(in) :: rs(*) double precision, intent(in) :: sorig(*) double precision, intent(in) :: r8heap(*) double precision, intent(in) :: r8heap1n(*) double precision, intent(inout) :: eigenvalue(*) double precision, intent(inout) :: c(*) call huzord_sch $(nbasis, c, rs, eigenvalue, r8heap, r8heap1n, iout, $nfroz, sorig, qmreg, 'schA', nal, qmreg_froz, .false. ) call huzord_sch $(nbasis, c(sqrsize+1), rs(sqrsize+1), eigenvalue(nbasis+1), $r8heap, r8heap1n,iout, nfroz, sorig, qmreg, 'schB', nal, $qmreg_froz, .false.) end subroutine sort_orbitals_for_sch_embedding ************************************************************************ ************************************************************************ subroutine huzord_sch(nbasis,c,rs,eig,r8heap2a_nbfXnbf, $r8heap2b_nbfXnbf,iout,nfroz,s,qmreg,route,nelec,qmreg_froz, $sort_virtuals) ************************************************************************ * Move frozen orbitals to the beginning in the case of Huzinaga eqs ************************************************************************ implicit none ! intent(in) character(len=4), intent(in) :: route character(len=8), intent(in) :: qmreg integer, intent(in) :: nelec integer, intent(in) :: nbasis integer, intent(in) :: nfroz integer, intent(in) :: iout integer, intent(in) :: qmreg_froz double precision, intent(in) :: s(nbasis,nbasis) double precision, intent(in) :: rs(nbasis,nbasis) double precision, intent(in), target :: r8heap2a_nbfXnbf $ (nbasis,nbasis) double precision, intent(in), target :: r8heap2b_nbfXnbf $ (nbasis,nbasis) logical, intent(in) :: sort_virtuals ! intent(inout) double precision, intent(inout) :: eig(nbasis) double precision, intent(inout) :: c(nbasis,nbasis) ! scr double precision, pointer :: work(:,:) =>null() double precision, pointer :: c_envXs(:,:) =>null() double precision, pointer :: MO_overlap(:) =>null() double precision, pointer :: rscr_nbf(:) =>null() ! external double precision, external :: ddot ! local integer :: fst_ind, lst_ind,i character(len=4) :: subsystem subsystem = route subsystem(1:3) = ' ' subsystem = trim(adjustl( subsystem ) ) call init_pointers call calc_MO_overlap if( subsystem .eq. 'A' ) then call sch_sort_MOs_subsystemA else if( subsystem .eq. 'B' ) then call sch_sort_MOs_subsystemB endif call check_for_intruder_states call print_warnings call set_nullpointers C contains subroutine init_pointers work => r8heap2a_nbfXnbf c_envXs => r8heap2b_nbfXnbf MO_overlap => r8heap2a_nbfXnbf(:,1) rscr_nbf => r8heap2a_nbfXnbf(:,2) work = 0.0d0 end subroutine init_pointers subroutine set_nullpointers work => null() c_envXs => null() MO_overlap => null() rscr_nbf => null() end subroutine set_nullpointers subroutine calc_MO_overlap integer :: iMO call dgemm('n','n',nbasis,nbasis,nbasis,1.d0,rs,nbasis,c, $nbasis,0.d0,work,nbasis) call dsymm('l','u',nbasis,nbasis,1.d0,s,nbasis,work, $nbasis,0.d0,c_envXs,nbasis) MO_overlap = 0.0d0 rscr_nbf = 0.0d0 do iMO = 1 , nbasis MO_overlap( iMO ) = ddot $ ( nbasis ,c_envXs(1,iMO) ,1 ,c(1,iMO),1) enddo end subroutine calc_MO_overlap subroutine sch_sort_MOs_subsystemA fst_ind = 1 lst_ind = nfroz call sort_orbitals $ ( get_largest_MO_overlap_ind,(/1,nbasis/)) call sort_orbitals $ ( get_highest_eig_ind, (/ 1 , lst_ind /) ) fst_ind = nfroz + 1 lst_ind = nbasis call sort_orbitals $ ( get_lowest_eig_ind, (/ 1 , lst_ind /) ) end subroutine sch_sort_MOs_subsystemA subroutine sch_sort_MOs_subsystemb fst_ind = 1 lst_ind = nfroz call sort_orbitals $ ( get_lowest_eig_ind, (/ 1 , lst_ind /) ) fst_ind = nfroz + 1 lst_ind = nelec call sort_orbitals $ ( get_largest_MO_overlap_ind,(/1,nbasis/)) call sort_orbitals $ ( get_highest_eig_ind, (/ 1 , lst_ind /) ) fst_ind = nelec + 1 lst_ind = nbasis call sort_orbitals $ ( get_lowest_eig_ind, (/ 1 , lst_ind /) ) end subroutine sch_sort_MOs_subsystemb subroutine sort_orbitals( func , func_range ) integer, external :: func integer, intent(in) :: func_range(2) integer :: i,ii double precision :: tmp i = fst_ind do while( i .le. lst_ind ) ii = func( i , func_range(2) ) if( i .eq. ii ) continue tmp=eig(i) eig(i)=eig(ii) eig(ii)=tmp rscr_nbf(1:nbasis)=c(1:nbasis,i) c(1:nbasis,i)=c(1:nbasis,ii) c(1:nbasis,ii)=rscr_nbf(1:nbasis) tmp=MO_overlap(i) MO_overlap(i)=MO_overlap(ii) MO_overlap(ii)=tmp i = i + 1 enddo end subroutine sort_orbitals pure integer function get_largest_MO_overlap_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = maxloc(MO_overlap(from:to),1)+from-1 end function get_largest_MO_overlap_ind pure integer function get_lowest_MO_overlap_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = minloc(MO_overlap(from:to),1)+from-1 end function get_lowest_MO_overlap_ind pure integer function get_highest_eig_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = maxloc(eig(from:to),1)+from-1 end function get_highest_eig_ind pure integer function get_lowest_eig_ind( from , to ) $ result( result_ind ) integer, intent(in) :: from integer, intent(in) :: to result_ind = minloc(eig(from:to),1)+from-1 end function get_lowest_eig_ind subroutine check_for_intruder_states integer :: nMOs integer :: maxlap_ind, minlap_ind double precision :: tmp, maxlap, minlap integer :: i fst_ind = merge( nfroz + 1 , 1 , subsystem.eq.'B') lst_ind = merge( nelec , nfroz , subsystem.eq.'B') minlap = $ MO_overlap( get_lowest_MO_overlap_ind( fst_ind, lst_ind ) ) fst_ind = merge( 1 , nfroz + 1 , subsystem.eq.'B') lst_ind = merge( nfroz , nelec , subsystem.eq.'B') maxlap_ind = get_largest_MO_overlap_ind( fst_ind, lst_ind ) maxlap = MO_overlap( maxlap_ind ) nMOs = merge( nbasis - nelec , nelec , sort_virtuals ) if( nfroz .ne. 0 .and. nfroz .ne. nMOs ) then if( maxlap .gt. 0.1d0 ) then write(iout,'(a)') $' Warning! Intruder states are detected!' write(iout,'(a)') $' This can cause convergence problems.' write(iout,'(a,f10.5,a,i8,a)') $' The largest overlap amongst intruder states: ', $maxlap,' (MO: ',maxlap_ind,')' endif tmp = maxlap - minlap if( abs(tmp) .lt. 0.1d0 ) then write(iout,*) write(iout,'(a)') $' Warning! Non-projector orbitals critically'// $' overlap with the projector.' write(iout,'(a)') $' This can cause convergence problems.' write(iout,'(a,f10.5)') $' Difference of orbital overlaps = ',abs(tmp) write(iout,*) endif endif end subroutine check_for_intruder_states subroutine print_warnings integer :: pMO fst_ind = merge( 1 , nfroz + 1 ,subsystem.eq.'A') lst_ind = merge( nfroz , nelec ,subsystem.eq.'A') do pMO = fst_ind , lst_ind if( eig( pMO ) .lt. 0.0d0 ) then call warning_msg( pMO ) endif enddo end subroutine print_warnings subroutine warning_msg( pMO ) integer, intent(in) :: pMO character(len=8) :: cscr write(cscr,'(i8)') pMO if( .not. sort_virtuals ) then write(iout,'(a,a,a)') $' Warning! The orbital energy of frozen occupied orbital ' $// trim(adjustl(cscr)) // ' is negative!' else write(iout,'(a)') $' Warning! The orbital energy of frozen virtual orbital ' $// trim(adjustl(cscr)) // ' is negative!' endif end subroutine warning_msg end subroutine huzord_sch ************************************************************************ subroutine save_f_embed(fock,symfock,r8heap,symtra,offset,lsa, $ nbasis,scrfile,itol) ************************************************************************ * Saving Ftilde[D^Atilde] for Huzinaga embedding gradient ************************************************************************ implicit none integer nbasis,scrfile,offset(*) double precision fock(*),symfock(*),r8heap(*),itol,symtra(*) logical lsa open(scrfile,file='FOCK_EMBED',form='unformatted') call mxto(symfock,symtra,r8heap,nbasis,fock,offset,lsa) call woeint(r8heap,r8heap,symfock,scrfile,itol,nbasis) close(scrfile) end subroutine save_f_embed ************************************************************************ ********************** ************************** ********************** UNUSED / UNTESTED ************************** ********************** ************************** ************************************************************************ subroutine fdmorb(nbasis,nb,no,natoms,sold,s,dold,d,fold,f,cold, $c,eig,scr1,scr2,scr3,minpfile,scrfile1,natrange,coord,verbosity, $cks,nocc) ************************************************************************ * Construct orbitals in the case of FDM embedding ************************************************************************ implicit none integer nbasis,nb,nv,no,natoms,ind(natoms),natrange(2,natoms),i integer minpfile,scrfile1,natomsn,verbosity,nsp,nocc real*8 sold(nbasis,nbasis),s(nb,nb),fold(nbasis,nbasis),f(nb,nb) real*8 dold(nbasis,nbasis),d(nb,nb),cold(nbasis,nbasis),c(nb,nb) real*8 scr1(nb,nb),scr2(nb,nb),scr3(nb,nb),eig(nb),coord(3,natoms) real*8 lshift,cks(nbasis,nbasis) parameter(lshift=1000000.d0) character*2 atsymbol(natoms) C Read atoms of embedded subsystem open(minpfile,file='MINP') call embedat(natoms,minpfile,ind,6,'embed ',5) close(minpfile) C Compress matrices call subsmat(natoms,nbasis,nb,fold,f,ind,natrange) call subsmat(natoms,nbasis,nb,dold,d,ind,natrange) call subsmat(natoms,nbasis,nb,sold,s,ind,natrange) C Construct PAOs call dsymm('l','l',nb,nb,-0.5d0,d,nb,s,nb,0.d0,c,nb) do i=1,nb c(i,i)=c(i,i)+1.d0 enddo C Remove linearly dependent PAOs and orthogonalize nv=nb-no call dsymm('l','l',nb,nb,1.d0,s,nb,c,nb,0.d0,scr1,nb) call dgemm('t','n',nb,nb,nb,1.d0,c,nb,scr1,nb,0.d0,scr2,nb) call dsyev('V','U',nb,scr2,nb,eig,scr3,3*nb**2,i) do i=no+1,nb call dscal(nb,1.d0/dsqrt(eig(i)),scr2(1,i),1) enddo call dgemm('n','n',nb,nv,nb,1.d0,c,nb,scr2(1,no+1),nb,0.d0, $scr1(1,no+1),nb) C Compress MO coefficients call subsvec(natoms,nbasis,nb,no,cold,scr1,ind,natrange) C Canonicalize orbitals call dsymm('l','u',nb,nb,1.d0,f,nb,scr1,nb,0.d0,scr3,nb) call dgemm('t','n',nb,nb,nb,1.d0,scr1,nb,scr3,nb,0.d0,scr2,nb) call dsyev('V','U',no,scr2,nb,eig,scr3,3*nb**2,i) call dsyev('V','U',nv,scr2(no+1,no+1),nb,eig(no+1),scr3,3*nb**2,i) call dgemm('n','n',nb,no,no,1.d0,scr1,nb,scr2,nb,0.d0,c,nb) call dgemm('n','n',nb,nv,nv,1.d0,scr1(1,no+1),nb,scr2(no+1,no+1), $nb,0.d0,c(1,no+1),nb) C Remove spurious occupied orbitals nsp=0 i=no+1 do while(i.le.nb.and.eig(i).lt.0.d0) i=i+1 nsp=nsp+1 enddo if(nsp.gt.0) then if(verbosity.ge.3) then write(6,*) write(6,*) 'Number of spurious occupied orbitals:',nsp write(6,*) $'Original orbital energies in embedded subsystem [au]:' write(6,"(7f10.6)") eig endif call dsyrk('u','n',nb,nsp,1.d0,c(1,no+1),nb,0.d0,scr1,nb) call dsymm('l','u',nb,nb,1.d0,scr1,nb,s,nb,0.d0,scr2,nb) call dsymm('l','u',nb,nb,lshift,s,nb,scr2,nb,1.d0,f,nb) call dcopy(nb*nsp,c(1,no+1),1,scr1,1) do i=1,nv-nsp c(1:nb,no+i)=c(1:nb,no+nsp+i) enddo call dcopy(nb*nsp,scr1,1,c(1,nb-nsp+1),1) call dcopy(nsp,eig(no+1),1,scr1,1) do i=1,nv-nsp eig(no+i)=eig(no+nsp+i) enddo call dcopy(nsp,scr1,1,eig(nb-nsp+1),1) do i=nb-nsp+1,nb eig(i)=eig(i)+lshift enddo endif C Save new MO coefficients and Fock matrix open(scrfile1,file='MOCOEF',form='unformatted') call wrtmo(scr1,scr1,c,scrfile1,0.d0,nb,nb) close(scrfile1) open(scrfile1,file='FOCK',form='unformatted') call woeint(scr1,scr1,f,scrfile1,0.d0,nb) write(scrfile1) (eig(i),i=1,nb) close(scrfile1) C Save coordinates call getvar('atsymbol ',atsymbol) call getvar('coord ',coord) natomsn=0 do i=1,natoms if(ind(i).gt.0) then natomsn=natomsn+1 atsymbol(natomsn)=atsymbol(i) coord(1:3,natomsn)=coord(1:3,i) endif enddo open(scrfile1,file='FDM',form='unformatted') write(scrfile1) natomsn,2*no write(scrfile1) atsymbol(1:natomsn) write(scrfile1) coord(1:3,1:natomsn) close(scrfile1) C if(verbosity.ge.3) then write(6,*) write(6,*) 'Orbital energies in embedded subsystem [au]:' write(6,"(7f10.6)") eig endif c szemet c call subsmat(natoms,nbasis,nb,sold,s,ind,natrange) c call dsymm('l','u',nb,nb,1.d0,s,nb,c,nb,0.d0,scr1,nb) c call dgemm('t','n',nb,nb,nb,1.d0,c,nb,scr1,nb,0.d0,d,nb) c write(6,*) 'd' c do i=1,nb c write(6,"(1000f10.6)") d(i,1:nb) c enddo c write(6,*) 'c' c do i=1,nb c write(6,"(1000f10.6)") c(i,1:nb) c enddo C return end C end module scf_embed