Merge pull request #1148 from edoapra/ses-pull

New SES cavity code
This commit is contained in:
NWChem: Open Source High-Performance Computational Chemistry 2025-09-01 12:02:18 -07:00 committed by GitHub
commit 2f7ab4a93a
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
22 changed files with 22715 additions and 313 deletions

View file

@ -736,6 +736,8 @@ let "myexit+=$?"
let "myexit+=$?"
./runtests.mpi.unix procs $np dft_nio_r2scan
let "myexit+=$?"
./runtests.mpi.unix procs $np cosmo_h2cco2_opt
let "myexit+=$?"
if [[ "$what" != "flaky" ]]; then
echo
echo "the number of failed tests is" $myexit

View file

@ -0,0 +1,81 @@
echo
start
geometry
o 0.00000000 -1.11909860 -0.09529177
o 0.00000000 1.11909860 -0.09529177
c 0.00000000 0.00000000 0.66838924
c 0.00000000 0.00000000 1.99717688
h 0.00000000 -0.92756722 2.53578926
h 0.00000000 0.92756722 2.53578926
end
charge -2
basis spherical
* library 6-31g*
end
title "new potential"
cosmo
potcorr true
ificos 1 ; minbem 3
do_gasphase f
end
driver
maxiter 99
clear
end
dft
xc hfexch
end
task dft optimize
title "old potential"
geometry
o 0.00000000 -1.11909860 -0.09529177
o 0.00000000 1.11909860 -0.09529177
c 0.00000000 0.00000000 0.66838924
c 0.00000000 0.00000000 1.99717688
h 0.00000000 -0.92756722 2.53578926
h 0.00000000 0.92756722 2.53578926
end
cosmo
potcorr false
do_gasphase f
ificos 1 ; minbem 3
print cosmo_energies
end
driver
clear
end
task scf optimize ignore
title "GEPOL cavity"
geometry
o 0.00000000 -1.11909860 -0.09529177
o 0.00000000 1.11909860 -0.09529177
c 0.00000000 0.00000000 0.66838924
c 0.00000000 0.00000000 1.99717688
h 0.00000000 -0.92756722 2.53578926
h 0.00000000 0.92756722 2.53578926
end
cosmo
do_gasphase f
cavity ses
ndiv xfine
end
unset cosmo:ptcor
driver
clear
end
task scf optimize ignore

File diff suppressed because it is too large Load diff

View file

@ -173,7 +173,8 @@ c
call ga_zero(g_dens)
call ga_copy(g_cdens,g_dens)
if (nopen.gt.0) call ga_add(1.d0,g_cdens,one,g_odens,g_dens)
call cosmo_charges_from_dmat(rtdb, basis, geom, ecosmo, odebug,
call cosmo_charges_from_dmat(rtdb, basis, geom,
& ecosmo, 0d0, 0d0, odebug,
& 1, ! 1 = closed shell
& g_dens, ! input density
& cosmo_file) ! cosmo charges file name

View file

@ -1152,7 +1152,8 @@ c check if we are doing dft (e.g. cphf case)
end do
call ga_copy(g_a_dens,g_dens(1)) ! alpha = 1
call ga_copy(g_b_dens,g_dens(2)) ! beta = 2
call cosmo_charges_from_dmat(rtdb, basis, geom, ecosmo, odebug,
call cosmo_charges_from_dmat(rtdb, basis, geom,
& ecosmo, 0d0, energy, odebug,
& 2, ! 2 = open shell
& g_dens, ! input density
& cosmo_file) ! cosmo charges file name

View file

@ -23,14 +23,11 @@ include ../config/makelib.h
ifeq ($(_FC),pgf90)
FOPTIONS += -Mnorecursive
endif
ifeq ($(FC),gfortran)
ifdef USE_OPENMP
FOPTIONS += -fno-recursive -fno-openmp
endif
endif
ifeq ($(_FC),gfortran)
ifdef USE_OPENMP
FOPTIONS += -fno-recursive -fno-openmp
ifneq ($(FLANG_NEW),true)
FOPTIONS += -fno-recursive -fno-openmp
endif
endif
endif
ifeq ($(FC),pathf90)

View file

@ -754,7 +754,7 @@ c if (oroot) write(luout,*) "calling cosmo_charges_from_dmat"
call cosmo_charges_from_dmat(rtdb,
& ihdl_bfao,
& ihdl_geom,
& ecosmo,
& ecosmo, 0d0, 0d0,
& odbug, ! debug
& 1, ! up and down are combined
& g_p_ao, ! input density

View file

@ -259,6 +259,7 @@ c == initialize DIIS call counter ==
c
gcds = 0.d0
egas = 0d0
esol = 0d0
do_cosmo_vem = 0
istep_cosmo_vem = 0
do_cosmo_smd = .false.
@ -1198,8 +1199,8 @@ c
cosmo_file = 'cosmo.xyz'
if (ipol.eq.2) ! separate components
& call ga_dadd(1.d0,g_dens(1),-1.d0,g_dens(2),g_dens(1))
call cosmo_charges_from_dmat(rtdb, ao_bas_han, geom,
& ecosmo, odbug, ipol, g_dens, cosmo_file)
call cosmo_charges_from_dmat(rtdb, ao_bas_han, geom, ecosmo,
& egas, etnew+enuc, odbug, ipol, g_dens, cosmo_file)
if (ipol.eq.2) ! reinstate total
& call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
end if ! cosmo check
@ -2192,8 +2193,8 @@ c
cosmo_file = 'cosmo.xyz'
if(ipol.eq.2) ! separate components
& call ga_dadd(1.d0,g_dens(1),-1.d0,g_dens(2),g_dens(1))
call cosmo_charges_from_dmat(rtdb, ao_bas_han, geom,
& ecosmo, odbug, ipol, g_dens, cosmo_file)
call cosmo_charges_from_dmat(rtdb, ao_bas_han, geom, ecosmo,
& egas, esol, odbug, ipol, g_dens, cosmo_file)
if (ipol.eq.2) ! reinstate total in g_dens(1)
& call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
end if ! cosmo_esp

View file

@ -274,7 +274,8 @@ c
call ga_zero(g_dens)
call ga_copy(g_cdens,g_dens)
if (nopen.gt.0) call ga_add(1.d0,g_cdens,one,g_odens,g_dens)
call cosmo_charges_from_dmat(rtdb, basis, geom, ecosmo, odebug,
call cosmo_charges_from_dmat(rtdb, basis, geom,
& ecosmo, 0d0, 0d0, odebug,
& 1, ! 1 = closed shell
& g_dens, ! input density
& cosmo_file) ! cosmo charges file name

View file

@ -348,7 +348,8 @@ c
end do
call ga_copy(g_a_dens,g_dens(1)) ! alpha = 1
call ga_copy(g_b_dens,g_dens(2)) ! beta = 2
call cosmo_charges_from_dmat(rtdb, basis, geom, ecosmo, odebug,
call cosmo_charges_from_dmat(rtdb, basis, geom,
& ecosmo, 0d0, energy, odebug,
& 2, ! 2 = open shell
& g_dens, ! input density
& cosmo_file) ! cosmo charges file name

View file

@ -1268,7 +1268,7 @@ c Calculate cosmo charges
if (cosmo_on.and.cosmo_phase.eq.2) then
cosmo_file = "cosmo.xyz"
call cosmo_charges_from_dmat(rtdb, ao_bas_han, geom,
& ecosmo, odbug, ipol, g_dens, cosmo_file)
& ecosmo, egas, etnew+enuc, odbug, ipol, g_dens, cosmo_file)
call cosmo_fock(ao_bas_han, g_fock)
call ga_dadd(one, g_dens(1), one, g_dens(2), g_dens(1))
Ecore = ga_ddot(g_dens(1), g_fock)

View file

@ -8,15 +8,18 @@
cosmo_smd_output.o\
cosmo_vem_output.o\
cosmo_lebedev.o\
cosmo_fibonacci.o
OBJ_OPTIMIZE = cosmo_input.o \
cosmo_fibonacci.o\
hnd_cosmo_lib.o hnd_coschg.o\
cosmo_input.o \
cosmo_initialize.o\
cosmo_tidy.o\
cosmo_def_radii.o\
cosmo_charges.o\
grad_hnd_cos.o hnd_cosmo_lib.o hnd_coschg.o\
gepol_lib.o\
grad_hnd_cos.o\
cosmo_thomson.o
OBJ_OPTIMIZE =
LIBRARY = libsolvation.a
@ -34,4 +37,7 @@ HEADERS = ../property/prop.fh ../property/atomdata.fh cosmoP.fh cosmo_params.fh
include ../config/makefile.h
include ../config/makelib.h
ifeq ($(FLANG_NEW),true)
FOPTIONS += -fno-automatic
endif
#FDEBUG += -O0

View file

@ -44,3 +44,12 @@ c
integer COSMO_LINEQ_ITERATIVE
parameter (COSMO_LINEQ_DIRECT = 0)
parameter (COSMO_LINEQ_ITERATIVE = 1)
c
c COSMO cavity type
c - COSMO_CAVITY_VDW : Use a van der Waals cavity
c - COSMO_CAVITY_SES : Use a solvent-excluded cavity
c
integer COSMO_CAVITY_VDW
integer COSMO_CAVITY_SES
parameter (COSMO_CAVITY_VDW = 0)
parameter (COSMO_CAVITY_SES = 1)

View file

@ -5,225 +5,11 @@ C>
C> \file cosmo_charges.F
C> The COSMO surface generator
C>
C> \brief Compute the COSMO charges
C>
subroutine cosmo_charges(rtdb,basis,geom,ecos,some,
& itype_wfn,vectors,nclosed,nopen,nbf,nmo)
c
implicit none
c
#include "errquit.fh"
#include "global.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "util.fh"
#include "mafdecls.fh"
#include "nwc_const.fh"
#include "stdio.fh"
#include "prop.fh"
#include "cosmoP.fh"
#include "cosmo_params.fh"
#include "util_params.fh"
c
integer rtdb !< [Input] The RTDB handle
integer basis !< [Input] The basis set handle
integer geom !< [Input] The geometry handle
integer nbf !< [Input] The number of basis functions
integer nmo !< [Input] The number of molecular orbitals
integer itype_wfn !< [Input] The wavefunction type
!< - 1: RHF type wavefunction
!< - 2: UHF type wavefunction
integer vectors(2) !< [Input] The GA handles for the SCF vectors
integer nclosed(2) !< [Input] The number of closed shell orbitals
!< (only used with RHF)
integer nopen(2) !< [Input] The number of open shell orbitals
c
double precision ecos !< [Output] The intra COSMO charges
!< interaction energy
logical some !< [Input] Do you want "some" additional output?
logical dbug, out, status
character*8 scfruntyp
c
character*255 cosmo_file
c
c
integer l_efcc, k_efcc, l_efcs, k_efcs, l_efcz, k_efcz
integer l_msrf, k_msrf
integer l_efciat, k_efciat
integer l_rad, k_rad, nrad
integer l_occ, k_occ
integer l_efczz, k_efczz
integer nefc, ief, i, nat
integer g_dens(3) ! 1: up, 2: down, 3: total
integer ndens,nvirt(2)
logical stat
c
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
c
dbug=.false..and.ga_nodeid().eq.0
out =.false..and.ga_nodeid().eq.0
out =out.or.dbug
some=some.or.out.and.ga_nodeid().eq.0
c
if(some) then
write(luout,9999)
endif
c
if(out) then
write(luout,*) 'in cosmo_charges ... geom = ',geom
endif
c
c ----- retrieve the number of atoms from -geom- -----
c
if (.not.geom_ncent(geom,nat))
& call errquit("cosmo_charges: geom_ncent failed",0,UERR)
c
c ----- read -efc- coordinates from -rtdb- -----
c
if(.not.rtdb_get(rtdb,'cosmo:nefc',mt_int,1 ,nefc))
& call errquit('cosmo_charges: rtdb get failed for nefc ',911,
& RTDB_ERR)
c
c ----- allocate memory for efc's -----
c
if(.not.ma_push_get(mt_dbl,nefc*3,'cosmo efcc',l_efcc,k_efcc))
& call errquit('cosmo_charges malloc k_efcc failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcs',l_efcs,k_efcs))
& call errquit('cosmo_charges malloc k_efcs failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'molsurf',l_msrf,k_msrf))
& call errquit('cosmo_charges malloc k_msrf failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcz',l_efcz,k_efcz))
& call errquit('cosmo_charges malloc k_efcz failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efczz',l_efczz,k_efczz))
& call errquit('cosmo_charges malloc k_efczz failed',911,MA_ERR)
if(.not.ma_push_get(mt_int,nefc,'cosmo efciat',l_efciat,k_efciat))
& call errquit('cosmo_charges malloc k_efciat failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nat,'cosmo rad',l_rad,k_rad))
& call errquit('cosmo_charges malloc k_rad failed',nat,MA_ERR)
c
if(.not.rtdb_get(rtdb,'cosmo:efcc',mt_dbl,3*nefc,dbl_mb(k_efcc)))
& call errquit('cosmo_charges: rtdb get failed for efcc ',912,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efcz',mt_dbl, nefc,dbl_mb(k_efcz)))
& call errquit('cosmo_charges: rtdb get failed for efcz ',913,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efcs',mt_dbl, nefc,dbl_mb(k_efcs)))
& call errquit('cosmo_charges: rtdb get failed for efcs ',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efczz',mt_dbl,nefc,dbl_mb(k_efczz)))
& call errquit('cosmo_charges: rtdb get failed for efczz',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efciat',
& mt_int,nefc,int_mb(k_efciat)))
& call errquit('cosmo_charges: rtdb get failed for efciat',914,
& RTDB_ERR)
call cosmo_def_radii(rtdb,geom,nat,dbl_mb(k_rad),.false.)
status = rtdb_get(rtdb,'cosmo:radius',mt_dbl, nat,dbl_mb(k_rad))
do i = 0, nat-1
dbl_mb(k_rad+i) = dbl_mb(k_rad+i)/cau2ang
enddo
c
if(out) then
write(luout,*) 'in cosmo_charges, nefc = ',nefc
do ief=1,nefc
write(luout,*) dbl_mb(k_efcc+(ief-1)*3),
& dbl_mb(k_efcc+(ief-1)*3+1),dbl_mb(k_efcc+(ief-1)*3+2)
enddo
endif
c
c ----- calculate the density matrices from molecular orbitals -----
cc
cc ----- get density matrix -----
cc
c scfruntyp='RHF'
c nocc=nclosed(1)+nopen(1)
c if (itype_wfn.eq.2) then
c scfruntyp='UHF'
c nocc=max(nocc,nclosed(2)+nopen(2))
c endif
c if(.not.ma_push_get(mt_dbl,nocc*2,'cosmo occ',l_occ,k_occ))
c & call errquit('cosmo_charges malloc k_occ failed',911,MA_ERR)
c do i=1,nocc*2
c dbl_mb(k_occ+i-1)=1.0d0
c enddo
c
c changes to make cosmo work for open shell DFT (MV)
c
if(.not.ma_push_get(mt_dbl,nbf*2,'cosmo occ',l_occ,k_occ))
& call errquit('cosmo_charges malloc k_occ failed',911,MA_ERR)
call dfill(2*nbf, 0.0d0, dbl_mb(k_occ), 1)
c
c Only need to set occupation numbers for UHF
c occupation numbers for RHF are done inside hnd_prop_dens_make
c
if (itype_wfn.eq.2) then
scfruntyp='UHF'
do i = 1, nopen(1)
dbl_mb(i-1+k_occ) = 1.0d0
enddo
do i = nbf+1, nbf+nopen(2)
dbl_mb(i-1+k_occ) = 1.0d0
enddo
else if (itype_wfn.eq.1) then
scfruntyp='RHF'
do i = 1, nclosed(1)
dbl_mb(i-1+k_occ) = 2.0d0
enddo
do i = nclosed(1)+1, nclosed(1)+nopen(1)
dbl_mb(i-1+k_occ) = 1.0d0
enddo
else
call errquit("unknown function type",0,0)
endif
c
c end of changes (MV)
c
call hnd_prop_dens_make(rtdb,geom,basis,nbf,nmo,nclosed,nopen,
& nvirt,scfruntyp,vectors,dbl_mb(k_occ),
& g_dens,ndens)
c
if(out) then
write(luout,*) 'in -cosmo_charges, wfntyp = ',itype_wfn
write(luout,*) 'in -cosmo_charges, scftyp = ',scfruntyp
if(dbug) then
write(luout,*) 'in -cosmo_charges, -da- ...'
call ga_print(g_dens(1))
if(scfruntyp.eq.'UHF ') then
write(luout,*) 'in -cosmo_charges, -db- ...'
call ga_print(g_dens(2))
endif
endif ! dbug
endif !out
c
c ----- get -cosmo- charges -----
c
call hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& dbl_mb(k_efcc),dbl_mb(k_efcs),dbl_mb(k_msrf),
& dbl_mb(k_efcz),dbl_mb(k_efczz),
& int_mb(k_efciat),dbl_mb(k_rad),ecos,cosmo_file)
c
c ----- release memory block -----
c
do i = 1, ndens
if (.not.ga_destroy(g_dens(i))) call
& errquit('cosmo_charges: ga_destroy failed g_dens',0,GA_ERR)
enddo
if(.not.ma_chop_stack(l_efcc))
& call errquit('cosmo_charges, ma_chop_stack of l_efcc failed',911,
& ma_err)
c
return
9999 format(/,10X,15(1H-),
1 /,10X,'-cosmo- charges',
2 /,10X,15(1H-))
end
c
C> \brief Compute the COSMO charges from an input density matrix
c
subroutine cosmo_charges_from_dmat(rtdb,basis,geom,ecos,some,
& ipol,g_dens_in,cosmo_file)
subroutine cosmo_charges_from_dmat(rtdb,basis,geom,
& ecos,egas,esol,
& some,ipol,g_dens_in,cosmo_file)
c
implicit none
c
@ -241,6 +27,8 @@ c
integer basis !< [Input] The basis set handle
integer geom !< [Input] The geometry handle
double precision ecos !< [Output] The intra COSMO charges, interaction energy
double precision egas !< [Input] The total gas energy
double precision esol !< [Input] The total solvated energy
logical some !< [Input] Do you want "some" additional output?
integer ipol !< [Input] 1: RHF type wavefunction, 2: UHF type wavefunction
integer g_dens_in(ipol) !< [Input] 1: up, 2: down (input density matrix to get the charges)
@ -374,7 +162,8 @@ c
call hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& dbl_mb(k_efcc),dbl_mb(k_efcs),dbl_mb(k_msrf),
& dbl_mb(k_efcz),dbl_mb(k_efczz),
& int_mb(k_efciat),dbl_mb(k_rad),ecos,cosmo_file)
& int_mb(k_efciat),dbl_mb(k_rad),
& ecos,egas,esol,cosmo_file)
c
c ----- release memory block -----
c

View file

@ -36,7 +36,7 @@ c
logical status
logical cosmo
c
integer iat, nrad, cosmo_screen
integer iat, nrad, cosmo_screen, cosmo_cavity
integer l_rad, k_rad, l_costag, k_costag, l_coszan, k_coszan
integer l_coscoor, k_coscoor
integer invnuc, nat, nefc, nmap
@ -251,9 +251,11 @@ c
if(.not.ma_pop_stack(l_costag)) call
& errquit('cosmo_init,ma_pop_stack of l_tag failed',911, MA_ERR)
c
do_cosmo_model = DO_COSMO_YK
if(.not.rtdb_get(rtdb,'cosmo:dielec',mt_dbl,1,dielec))
D dielec =78.4d0 ! default water dielectric
cosmo_cavity = COSMO_CAVITY_VDW
if(.not.rtdb_get(rtdb,'cosmo:cavity',mt_int,1,cosmo_cavity))
D cosmo_cavity = COSMO_CAVITY_VDW
cosmo_screen = COSMO_SCREEN_ST
cosmo_sccor = COSMO_SCCOR_LAGRA ! surface charge correction adapted for geometry optimzations
minbem =2
@ -285,8 +287,9 @@ c
c
status = rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo)
status = rtdb_get(rtdb,'cosmo:rsolv',mt_dbl,1,rsolv)
status = rtdb_get(rtdb,'cosmo:do_cosmo_model',mt_int,1,
+ do_cosmo_model)
if(.not.rtdb_get(rtdb,'cosmo:do_cosmo_model',mt_int,1,
+ do_cosmo_model))
C do_cosmo_model=DO_COSMO_YK
if (do_cosmo_model.eq.DO_COSMO_KS) then
cosmo_sccor = COSMO_SCCOR_SCALE
cosmo_screen = COSMO_SCREEN_KS
@ -296,6 +299,8 @@ c
endif
status = rtdb_get(rtdb,'cosmo:screen',mt_int,1,cosmo_screen)
status = rtdb_get(rtdb,'cosmo:sccor',mt_int,1,cosmo_sccor)
if(.not.rtdb_get(rtdb,'cosmo:ptcor',mt_log,1,cosmo_ptcor))
I cosmo_ptcor = .false. ! correct potentials in case of outlying charge corrections
status = rtdb_get(rtdb,'cosmo:minbem',mt_int,1,minbem)
status = rtdb_get(rtdb,'cosmo:maxbem',mt_int,1,maxbem)
status = rtdb_get(rtdb,'cosmo:ificos',mt_int,1,ificos)
@ -506,6 +511,7 @@ c
write(Luout,9999)
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(Luout,9970)
if (cosmo_cavity.eq.COSMO_CAVITY_SES) write (Luout,8870)
else if (do_cosmo_model.eq.DO_COSMO_YK) then
write(Luout,9971)
endif
@ -608,15 +614,17 @@ c
if(osome.and.me.eq.0) then
write(luout,9993) lineq
if ((iangleb.lt.1).and.(fibonacci.lt.1)) then
write(luout,9991) minbem
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(luout,9992) maxbem
endif
if(ificos.eq.0) then
write(luout,9990)
else
write(luout,9989)
endif
if (cosmo_cavity.ne.COSMO_CAVITY_SES) then
write(luout,9991) minbem
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(luout,9992) maxbem
endif
if(ificos.eq.0) then
write(luout,9990)
else
write(luout,9989)
endif
end if
elseif (fibonacci.lt.1) then
write(luout,9969)
write(luout,9968) iangleb
@ -750,6 +758,8 @@ c
9972 format(' surface charge correction = ',a)
9971 format(' Cosmo: York-Karplus, doi: 10.1021/jp992097l')
9970 format(' Cosmo: Klamt-Schuurmann, doi: 10.1039/P29930000799')
c TODO: describe the GEPOL modifications and publish a new article
8870 format(' SES surface, doi: 10.1002/jcc.540151009, ??')
9969 format(' Using -Lebedev- grid')
9968 format(' -iangleb',I4)

View file

@ -72,6 +72,7 @@ c
integer iat, nfield, nrad, irad
integer iscren, ifscrn, minbem, maxbem, ificos, lineq
integer cosmo_screen, cosmo_sccor
logical cosmo_ptcor
integer iangleb
integer fibonacci
integer thomson
@ -79,6 +80,9 @@ c
double precision dielec,dielecinf,rsolv,rad,zeta,gammas,swtol
character*255 field
logical do_cosmo
integer cosmo_cavity
integer cosmo_ndiv
double precision mfactor
logical do_cosmo_file
logical do_gasphase
logical status
@ -136,6 +140,8 @@ c
solvname_long = 'water'
cosmo_screen = COSMO_SCREEN_ST
cosmo_sccor = COSMO_SCCOR_LAGRA ! surface charge correction adapted for geometry optimzations
cosmo_ndiv = 2 ! used with the SES cavity
mfactor =1.00d+00 ! used with the SES cavity
minbem =2
maxbem =4 ! used with the Klamt-Schuurmann model
ificos =0
@ -311,21 +317,25 @@ c
else if(inp_compare(.false.,'off',field)) then
do_cosmo=.false.
c
c --- do_cosmo_ks: switch to Klamt-Schuurmann model
c --- do_cosmo_ks: switch to Klamt-Schuurmann model; also makes cavity=ses
c
else if(inp_compare(.false.,'do_cosmo_ks',field)) then
do_cosmo_model=DO_COSMO_KS
if (.not. rtdb_put(rtdb,'cosmo:do_cosmo_model',mt_int,1,
$ do_cosmo_model))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
cosmo_cavity=COSMO_CAVITY_SES
if (.not.(rtdb_put(rtdb,'cosmo:do_cosmo_model',mt_int,1,
$ do_cosmo_model)
$ .and.rtdb_put(rtdb,'cosmo:cavity',mt_int,1,cosmo_cavity))
$ ) call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- do_cosmo_yk: switch to York-Karplus model
c --- do_cosmo_yk: switch to York-Karplus model; also makes cavity=vdw
c
else if(inp_compare(.false.,'do_cosmo_yk',field)) then
do_cosmo_model=DO_COSMO_YK
if (.not. rtdb_put(rtdb,'cosmo:do_cosmo_model',mt_int,1,
$ do_cosmo_model))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
cosmo_cavity=COSMO_CAVITY_VDW
if (.not.(rtdb_put(rtdb,'cosmo:do_cosmo_model',mt_int,1,
$ do_cosmo_model)
$ .and.rtdb_put(rtdb,'cosmo:cavity',mt_int,1,cosmo_cavity)
$ )) call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- do_cosmo_file flag
c
@ -333,8 +343,16 @@ c
do_cosmo_file=.true.
if (.not.
$ rtdb_put(rtdb,'cosmo:do_cosmo_file',mt_log,1,do_cosmo_file))
$ call errquit('cosmo_input: cannot put do_cosmo_file in rtdb',
$ 0,rtdb_err)
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- info: Information line to be printed in the cosmo file as a comment.
c Maximum lenth is 255 characters.
c
else if(inp_compare(.false.,'info',field)) then
if (.not.inp_a(token))
$ call errquit('cosmo_input: info',911, INPUT_ERR)
if (.not.(rtdb_put(rtdb,'cosmo:info',mt_byte,255,token)))
$ call errquit('cosmo_input: rtdb put failed',0,RTDB_ERR)
c
c --- rsolv: Used in Klamt-Schuurmann, but not in York-Karplus
c --- approach
@ -345,6 +363,84 @@ c
if (.not. rtdb_put(rtdb,'cosmo:rsolv',mt_dbl,1,rsolv))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- CAVITY: cavity construction method
c CAVITY VDW : van der Waals cavity
c CAVITY SES : solvent-excluded-cavity
c
c default for each model:
c - CAVITY SES : with the Klamt-Schuurmann model
c - CAVITY VDW : with the York-Karplus model
c
else if(inp_compare(.false.,'cavity',field)) then
if (.not.inp_a(token))
$ call errquit('cosmo_input: cavity',911, INPUT_ERR)
if (inp_compare(.false.,'vdw',token)) then
cosmo_cavity = COSMO_CAVITY_VDW
else if (inp_compare(.false.,'ses',token)) then
cosmo_cavity = COSMO_CAVITY_SES
c use KS model by default when using SES cavity
do_cosmo_model=DO_COSMO_KS
c correct potential when using SES
cosmo_ptcor=.true.
else
write(luout,*)"invalid option for CAVITY, valid options are:"
write(luout,*)"- VDW for a van der Waals cavity"
write(luout,*)"- SES for a solvent-excluded-cavity"
call util_flush(luout)
call errquit("cosmo_input: invalid option for CAVITY",911,
$ INPUT_ERR)
endif
if (.not. rtdb_put(rtdb,'cosmo:cavity',mt_int,1,cosmo_cavity))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
if (.not. rtdb_put(rtdb,'cosmo:do_cosmo_model',mt_int,1,
C do_cosmo_model))
$ call errquit('cosmo_input: rtdb put failed',912, RTDB_ERR)
if (.not.
$ rtdb_put(rtdb,'cosmo:ptcor',mt_log,1,cosmo_ptcor))
$ call errquit('cosmo_input: rtdb put failed',913, RTDB_ERR)
c
c --- ndiv: specifies the division level for the triangles on
c --- sphere surface when a SES cavity is used.
c Valid values are:
c - coarse: 60 triangles per sphere
c - fine : 240 (default)
c - xfine : 960
c - huge : 3840
c - xxl :15360
c
else if(inp_compare(.false.,'ndiv',field)) then
if (.not.inp_a(token))
$ call errquit('cosmo_input: ndiv',911, INPUT_ERR)
if (inp_compare(.false.,'coarse',token)) then
cosmo_ndiv = 1
else if (inp_compare(.false.,'fine',token)) then
cosmo_ndiv = 2
else if (inp_compare(.false.,'xfine',token)) then
cosmo_ndiv = 3
else if (inp_compare(.false.,'huge',token)) then
cosmo_ndiv = 4
else if (inp_compare(.false.,'xxl',token)) then
cosmo_ndiv = 5
else
write(luout,*)"invalid option for NDIV, valid options are:"
write(luout,*)" coarse, fine (default), xfine"
call util_flush(luout)
call errquit("cosmo_input: invalid option for NDIV",911,
$ INPUT_ERR)
endif
if (.not. rtdb_put(rtdb,'cosmo:ndiv',mt_int,1,cosmo_ndiv))
$ call errquit('cosmo_input: rtdb put failed',911,RTDB_ERR)
c
c --- mfactor: Merging factor for segments that are too close.
c Any two segments closer than f(rsolv,ndiv) * mfactor are merged into
c a single segment. The default value is 1.0.
c
else if(inp_compare(.false.,'mfactor',field)) then
if(.not. inp_f(mfactor))
$ call errquit('cosmo_input: mfactor',911, INPUT_ERR)
if (.not. rtdb_put(rtdb,'cosmo:mfactor',mt_dbl,1,mfactor))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- parameter file
c
else if(inp_compare(.false.,'parameters',field)) then
@ -481,6 +577,15 @@ c
if (.not. rtdb_put(rtdb,'cosmo:sccor',mt_int,1,cosmo_sccor))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- potcorr flag: use the potential
c correction associated with charge corrections
c
else if(inp_compare(.false.,'potcorr',field)) then
if (.not. inp_l(cosmo_ptcor)) cosmo_ptcor=.true.
if (.not.
$ rtdb_put(rtdb,'cosmo:ptcor',mt_log,1,cosmo_ptcor))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- minbem
c
else if(inp_compare(.false.,'minbem',field)) then

View file

@ -5,9 +5,10 @@ c
double precision adiag,dsurf,dvol,srfmol,volmol,ptspatm
integer do_cosmo_model,lineq,minbem,maxbem,ificos,ifscrn,
& cosmo_sccor,iangleb,fibonacci,thomson
logical cosmo_ptcor
common/hnd_cospar/dielec,dielecinf,screen,rsolv,zeta,gammas,swtol
common/hnd_cosmod/do_cosmo_model,lineq,minbem,maxbem,
I ificos,ifscrn,
& cosmo_sccor,iangleb,fibonacci,thomson
& cosmo_ptcor,cosmo_sccor,iangleb,fibonacci,thomson
common/hnd_cosdat/adiag,dsurf,dvol,srfmol,volmol,ptspatm

2593
src/solvation/gepol_lib.F Normal file

File diff suppressed because it is too large Load diff

View file

@ -445,10 +445,22 @@ C> "Envelope theorems for arbitrary choice sets",
C> <i>Econometrica</i> (2002) <b>70</b>, 583-601, DOI:
C> <a href="https://doi.org/10.1111/1468-0262.00296">
C> 10.1111/1468-0262.00296</a>.
C>
C> [10] A. Klamt, V. Jonas, T. Burger, J. C. W. Lohrenz,
C> "Refinement and Parametrization of COSMO-RS",
C> <i>J. Phys. Chem. A</i> (1998) <b>102</b>, 1089-5639, DOI:
C> <a href="https://doi.org/10.1021/jp980017s">
C> 10.1021/jp980017s</a>.
C>
C> [11] A. Klamt, V. Jonas,
C> "Treatment of the outlying charge in continuum solvation models",
C> <i>J. Chem. Phys.</i> (1996) <b>105</b>, pp 99729981, DOI:
C> <a href="https://doi.org/10.1063/1.472829">
C> 10.1063/1.472829</a>.
C>
subroutine hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& efcc,efcs,msrf,efcz,efczz,efciat,ratm,
& ecos,cosmo_file_in)
& ecos,egas,esol,cosmo_file_in)
implicit none
#include "errquit.fh"
#include "global.fh"
@ -489,11 +501,15 @@ c
double precision ratm(nat) !< [Input] the atom radii
double precision ecos !< [Output] the energy contribution due to
!< the COSMO charges
double precision egas !< [Input] the total gas energy
double precision esol !< [Input] the total solvated energy
logical status
logical do_cosmo_file
integer cosmo_cavity
c
logical dbug,more,out,direct,noall,all,elec,nucl,iefc_done
character*16 at_tag
character*255 cosmo_info
integer istrlen
character*255 cosmo_file_in
character*(nw_max_path_len) cosmo_file
@ -528,13 +544,13 @@ c integer cosmo_sccor ! do correction?
integer nodcmp ! flag specifying how to handle errors
integer need ! the amount of memory needed
c
integer l_epot, l_xyzpt, l_zanpt ! memory handles
integer k_epot, k_xyzpt, k_zanpt ! memory offsets
integer l_epot, l_xyzpt, l_zanpt, l_costag ! memory handles
integer k_epot, k_xyzpt, k_zanpt, k_costag ! memory offsets
c
double precision charge ! the total QM region charge
double precision chgnuc ! the total nuclear charge
integer nelec ! the total number of electrons
double precision chgfac ! scale factor for COSMO charges
double precision chgfac, chgfacn, chgface ! scale factor for COSMO charges
double precision chgcos ! the total COSMO surface charge
double precision chgcvg ! the convergence of the COSMO charges
double precision chgina ! the inv(A) COSMO charge
@ -556,7 +572,6 @@ c
double precision solnrg ! solvation energy
double precision dlambda ! lambda (surface charge correction)
double precision elambda ! lambda dependent energy term
double precision ecos2
double precision pi
double precision zero, pt5, one, two ! constants
data zero /0.0d+00/
@ -719,9 +734,11 @@ c
& call errquit('hnd_coschg: ma failed',911,MA_ERR)
if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt))
& call errquit('hnd_coschg: ma failed',911,MA_ERR)
if (.not. ma_push_get(mt_byte,nat*16,'tags',l_costag,k_costag))
& call errquit('hnd_coschg: k_costag failed',911,MA_ERR)
do iat=1,nat
if(.not.geom_cent_get(geom,iat,at_tag,dbl_mb(k_xyzpt+3*(iat-1)),
& dbl_mb(k_zanpt+iat-1))) call
if(.not.geom_cent_get(geom,iat,byte_mb(k_costag+16*(iat-1)),
& dbl_mb(k_xyzpt+3*(iat-1)),dbl_mb(k_zanpt+iat-1))) call
& errquit('hnd_coschg: geom_cent_get',911,GEOM_ERR)
enddo ! iat
c
@ -1005,9 +1022,11 @@ c
else if (nucl) then
errcos=chgnuc-chgcos
chgfac=chgnuc/chgcos
chgfacn = chgfac
else if (elec) then
errcos=-dble(nelec)-chgcos
chgfac=-dble(nelec)/chgcos
chgface = chgfac
endif
call dscal(nefc,chgfac,dbl_mb(ieq),1)
else if (cosmo_sccor.eq.COSMO_SCCOR_LAGRA) then
@ -1042,6 +1061,23 @@ c
enddo
endif
c
c ----- update potentials due to outlying charge correction
c
if (cosmo_sccor.ne.COSMO_SCCOR_NO.and.cosmo_ptcor) then
c If outlying charge corrections have been applied, the total
c potential should be updated, pot = -A*q, see Eq. 12 of Ref. [11].
c
c The total potential is stored at i10, total point charges at i20.
call hnd_cosaxd(nat,dbl_mb(i20),dbl_mb(i10),nefc,efcc,
+ efcs,efczz,efciat,ratm)
c The nuclear potential should not be changed by the outlying
c charges, so the elec potential should be adjusted.
c Elec potential (i12) = total (i10) - nuc (i11)
do ief=1,nefc
dbl_mb(ief+i12-1) = dbl_mb(ief+i10-1) - dbl_mb(ief+i11-1)
enddo
endif
c
c ----- charge screening due to the dielectric medium -----
c
c ----- set screening factor -----
@ -1106,6 +1142,13 @@ c surface charge corrections are used the results should be
c identical. If surface charge corrections are employed using
c scaling factors or Lagrangian multipliers then additional
c correction terms are needed to ensure the expressions agree.
c
c A consistent correction is to recalculate the total potential
c with the corrected charges (see Ref. [11] and cosmo_ptcor above).
c After that, adjust the elec potential accordingly, with a
c constant nuclear potential. This makes both formulations
c agree and also makes scale and Lagrangian corrections
c to produce similar results, elambda must be disregarded.
c
call hnd_cos_energy(nat,nefc,chgfac,efcc,efcs,efcz,efczz,efciat,
& ratm,dbl_mb(k_xyzpt),dbl_mb(k_zanpt),
@ -1131,22 +1174,22 @@ c
allefc = ddot(nefc,efcz,1,dbl_mb(i10),1)
atmefc = ddot(nefc,efcz,1,dbl_mb(i11),1)
elcefc = ddot(nefc,efcz,1,dbl_mb(i12),1)
solnrg= pt5* allefc + elambda
c solnrg= pt5* allefc + elambda
solnrg = pt5*allefc
ecos = pt5*(atmefc-elcefc)
c COSMO contribution Alternative 2 (Ref.[1] Eq.(11))
c ecos = pt5*(atmefc-elcefc)+elambda
if (do_cosmo_smd) then
ecos = pt5*(atmefc-elcefc)
else
c ecos = efcefc + atmefc
ecos = pt5*(atmefc-elcefc) + elambda
endif
c add lagrangian multiplier correction to the energy
c when the potential is not corrected
c
if(.not.cosmo_ptcor) ecos=ecos+elambda
if (oprint_energies) then
write(luout,*)'Alternative 2'
write(luout,9991) atmefc
write(luout,9990) elcefc
write(luout,9989) allefc
write(luout,9988) solnrg
write(luout,9989) allefc
write(luout,9994) ecos
endif
@ -1223,6 +1266,8 @@ c
if (.not. rtdb_get(rtdb,'cosmo:do_cosmo_file',mt_log,1,
$ do_cosmo_file)) do_cosmo_file=.false.
if (do_cosmo_file) then
if (.not. rtdb_get(rtdb,'cosmo:cavity',mt_int,1,
$ cosmo_cavity)) cosmo_cavity=COSMO_CAVITY_VDW
cosmo_file_in = "cosmo"
call util_file_name(cosmo_file_in,.false.,.false.,cosmo_file)
call util_file_name_resolve(cosmo_file,.false.)
@ -1235,37 +1280,78 @@ c
write(*,*) "printing cosmo file for COSMO-RS/SAC",
+ cosmo_file
end if
if (rtdb_get(rtdb,'cosmo:info',mt_byte,255,cosmo_info))
+ write(fn,"(a)") '# ' // trim(cosmo_info)
write(fn,"(a)") "$cosmo"
write(fn,*) " epsilon= ", dielec
write(fn,*) " rsolv= ", rsolv
write(fn,"(a,f6.3)") " screen= ", screen
if(screen.eq.one) then
write(fn,*) " epsilon= infinity"
else
write(fn,"(a,f6.3)") " epsilon= ", dielec
endif
if(cosmo_cavity.eq.COSMO_CAVITY_SES) then
write(fn,*) " cavity= SES"
write(fn,"(a,f6.3)") " rsolv= ", rsolv
else
write(fn,*) " cavity= vdW"
endif
if(do_cosmo_model.eq.DO_COSMO_YK) then
write(fn,*) " model= YK"
write(fn,"(a,f6.3)") " zeta= ", zeta
else
write(fn,*) " model= KS"
endif
if (cosmo_sccor.eq.COSMO_SCCOR_NO) then
write(fn,*) " charge_correction= no"
else if (cosmo_sccor.eq.COSMO_SCCOR_LAGRA) then
write(fn,*) " charge_correction= lagrangian"
else if (cosmo_sccor.eq.COSMO_SCCOR_SCALE) then
write(fn,*) " charge_correction= scale"
write(fn,"(a,f12.9)") " charge_corr_elec= ", chgface
write(fn,"(a,f12.9)") " charge_corr_nucl= ", chgfacn
end if
c
write(fn,"(a)") "$cosmo_data"
write(fn,*) " nps= ", nefc
write(fn,*) " area= ", srfmol/(cau2ang**2)
write(fn,*) " volume= ", volmol/(cau2ang**3)
write(fn,"(a,i6)") " nps= ", nefc
write(fn,"(a,f10.2)") " area= ", srfmol/(cau2ang**2)
write(fn,"(a,f10.2)") " volume= ", volmol/(cau2ang**3)
c
write(fn,"(a)") "$coord_rad"
write(fn,"(a)") "#atom x y z [a.u.] element radius [A]"
do i = 1,nat
write(fn,*) i,
call util_set_ma_char_arrays(16,
+ byte_mb(k_costag+16*(i-1)),at_tag)
write(fn,"(i6,3f15.9,2x,a,f6.3)") i,
+ dbl_mb(k_xyzpt +3*(i-1)),
+ dbl_mb(k_xyzpt+1+3*(i-1)),
+ dbl_mb(k_xyzpt+2+3*(i-1)),
+ int(dbl_mb(k_zanpt+i-1)),
+ at_tag,
+ ratm(i)*cau2ang
enddo ! i
c
c
write(fn,"(a)") "$energy"
if(dabs(egas).gt.0d0)
+ write(fn,"(a,f18.10)") " total_gas= ", egas
if(dabs(esol).gt.0d0)
+ write(fn,"(a,f18.10)") " total_solvated= ", esol
c
c Dielectric energy = 0.5 * sum(potential_i * q_i), Eq. (5) of Ref. [10].
c This energy is expected to correlate well with the solvation energy.
write(fn,"(a,f18.10)") " diel= ", pt5*allefc
c
write(fn,"(a)") "$segment_information"
write(fn,"(a)") "# n atom x y z charge area charge/area"
write(fn,"(a)") "# coordinates [a.u.], area [A2]"
write(fn,"(a)")
+ "# n atom x y z charge area charge/area potential"
write(fn,"(a)")
+ "# coordinates and potential [a.u.], area [A2]"
do ief=1,nefc
efcsi = efcs(ief)*(cau2ang**2)
if (do_cosmo_model.eq.DO_COSMO_YK)
+ efcsi = msrf(ief)*(cau2ang**2)
write(fn,*) ief, efciat(ief),
write(fn,"(2i6,7f15.9)") ief, efciat(ief),
+ efcc(1,ief), efcc(2,ief), efcc(3,ief),
+ efcz(ief), efcsi, efcz(ief)/efcsi
+ efcz(ief), efcsi, efcz(ief)/efcsi,
+ dbl_mb(i10+ief-1)
end do
close(fn)
end if

View file

@ -190,6 +190,8 @@ c
double precision wleb(mxfac)
double precision energy
c
integer cosmo_cavity
double precision mfactor
logical oprint_out
logical oprint_more
logical oprint_debug
@ -204,6 +206,11 @@ c
D (ga_nodeid().eq.0)
oprint_more = util_print("cosmo_cossrf_more",print_never).and.
D (ga_nodeid().eq.0)
c
if (.not. rtdb_get(rtdb,'cosmo:cavity',mt_int,1,
$ cosmo_cavity)) cosmo_cavity=COSMO_CAVITY_VDW
if (.not. rtdb_get(rtdb,'cosmo:mfactor',mt_dbl,1,
$ mfactor)) mfactor=1.d0
c
c ----- approximate sphere with segments and points -----
c
@ -267,7 +274,8 @@ c
if(radius(iat).eq.0.0d0) then
ratm(iat)=0.0d0
else
if (do_cosmo_model.eq.DO_COSMO_KS) then
if (do_cosmo_model.eq.DO_COSMO_KS.and.
+ cosmo_cavity.ne.COSMO_CAVITY_SES) then
ratm(iat)=(radius(iat)+rsolv)/cau2ang
else
ratm(iat)=radius(iat)/cau2ang
@ -280,8 +288,8 @@ c
call hnd_cossas(nat,xyzatm,ratm,mxatm,
1 nspa,nppa,xyzspa,ijkspa,
2 nseg,nfac,xyzseg,ijkseg,insseg,
3 numpps,xyzff,wleb,mxfac,geom,rtdb)
2 nseg,nfac,xyzseg,ijkseg,insseg,numpps,
3 xyzff,wleb,mxfac,geom,cosmo_cavity,mfactor,rtdb)
c
return
9999 format(' nseg,nfac,ndiv=nfac/nseg,dsurf,dvol = ',3i7,2f10.6)
@ -373,9 +381,9 @@ C> 10.1021/jp992097l</a>.
C>
subroutine hnd_cossas(nat,xyzatm,ratm,mxatom,
1 nspa,nppa,xyzspa,ijkspa,
2 nseg,nfac,xyzseg,ijkseg,insseg,
3 numpps,xyzff,wleb,mxface,geom,
4 rtdb)
2 nseg,nfac,xyzseg,ijkseg,insseg,numpps,
3 xyzff,wleb,mxface,geom,cosmo_cavity,
4 mfactor,rtdb)
implicit none
#include "cosmo_params.fh"
#include "errquit.fh"
@ -424,15 +432,24 @@ c
integer numpps( mxface,mxatom)
double precision xyzff( mxface,mxatom)
double precision wleb( mxface)
double precision zero, one
data one /1.0d+00/
double precision zero, one, two
data zero/0.0d00/,one/1.0d00/,two/2.0d00/
integer l_efcc, k_efcc, l_efcs, k_efcs, l_efcz, k_efcz
integer l_efclb, k_efclb, k_efciat, l_efciat
integer l_efczz, k_efczz
integer l_msrf,k_msrf
c SES variables
integer cosmo_cavity
double precision mfactor, asum, dmin
real*4 ofac,rd,dvec,rmim
INTEGER*4 ndiv, gseg
double precision segarea(mxface)
integer segatom(mxface)
double precision garea, gvolume
c
double precision ratm_real,dij,dum,cavdsp,pi,zetai,zetaii
integer m,mfac,mseg
integer nefc,iat,jat,npp,i,iseg,ifac,ief,ipp
integer nefc,iat,jat,npp,i,iseg,jseg,ifac,ief,ipp
integer l_dum1, l_dum2, l_attag, k_dum1, k_dum2, k_attag
character(len=16) aname
@ -470,12 +487,22 @@ c
c
oprint_debug = util_print("cosmo_cossas_debug",print_debug).and.
D (ga_nodeid().eq.0)
dbug=oprint_debug
oprint_molsurf = util_print("cosmo_mol_surface",print_never).and.
D (ga_nodeid().eq.0)
pi = acos(-1.0d0)
c
if(ga_nodeid().eq.0) then
write(luout,9999)
if(cosmo_cavity.eq.COSMO_CAVITY_SES) then
write(luout,8888)
else
c YK model is currently using a vdW surface and not a SAS
c as described in the original paper.
c The KS model with a vdW cavity has a defect: it produces an
c incomplete vdW surface, leaving solvent-excluded areas
c open, small values of rsolv are needed.
write(luout,9999)
endif
endif
c
c ----- print atomic centers -----
@ -483,11 +510,12 @@ c
if(ga_nodeid().eq.0) then
write(luout,9998)
do iat=1,nat
if (do_cosmo_model.eq.DO_COSMO_KS) then
if (do_cosmo_model.eq.DO_COSMO_KS.and.
+ cosmo_cavity.ne.COSMO_CAVITY_SES) then
write(luout,9997) iat,xyzatm(1,iat),xyzatm(2,iat),
1 xyzatm(3,iat),
2 (ratm(iat)*cau2ang-rsolv)
else if (do_cosmo_model.eq.DO_COSMO_YK) then
else
write(luout,9997) iat,xyzatm(1,iat),xyzatm(2,iat),
1 xyzatm(3,iat),
2 (ratm(iat)*cau2ang)
@ -669,7 +697,7 @@ c
c
enddo
c
if(ga_nodeid().eq.0) then
if(ga_nodeid().eq.0.and.cosmo_cavity.ne.COSMO_CAVITY_SES) then
write(luout,9985) nseg,nfac
write(luout,9992)
do iat=1,nat
@ -680,7 +708,7 @@ c
write(luout,9991) iat,nspa(iat),nppa(iat),npp
enddo
endif
if(oprint_debug) then
if(oprint_debug.and.cosmo_cavity.ne.COSMO_CAVITY_SES) then
write(luout,9987)
do iat=1,nat
do iseg=1,nseg
@ -701,6 +729,25 @@ c
endif
enddo
c
c Construct a precise SES surface, using a modified GEPOL93.
c
if (cosmo_cavity.eq.COSMO_CAVITY_SES) then
c We are sending all arguments in A.U.
rd = rsolv/cau2ang
rmim = rd/2
ofac = 0.9 ! to reduce the concave regions
dvec = 1.0 ! segment vector direction, not relevant
if(.not.rtdb_get(rtdb,'cosmo:ndiv',mt_int,1,ndiv)) ndiv=2
call gepol('ESURF',.false.,.false.,.false.,.false.,'TIT',0,
& rmim,'','','',ofac,rd,ndiv,dvec,
& .true.,.false.,
& nat,ratm,xyzatm,mxface,
& gseg,xyzseg,segarea,segatom,srfmol,volmol,.false.)
nefc = gseg
nseg = gseg
end if
c
c
c Allocate memory for point charges
c
if(.not.ma_push_get(mt_dbl,nefc*3,'cosmo efcc',l_efcc,k_efcc))
@ -736,6 +783,72 @@ c ----- save coordinates of surface points -----
c save segment surfaces
c save segment to atom mapping
c
if(cosmo_cavity.eq.COSMO_CAVITY_SES) then
c Save gepol results, but keep the segments sorted by atom
c so grad_hdn_cos can work properly.
c
c The minimum distance we will accept between 2 segments,
c this is slight lower than the minimum distance for a rsolv sphere.
dmin = rsolv/cau2ang * 0.3 * mfactor
if (ndiv.eq.2) dmin = dmin/2.0
if (ndiv.ge.3) dmin = dmin/4.0
c
c Checking for overlapping segments to merge.
if (do_cosmo_model.ne.DO_COSMO_YK.and.
A dmin.gt.zero) then
iseg = 0
do while (iseg <= nseg)
iseg = iseg + 1
if(segarea(iseg).le.zero) cycle
do jseg=iseg+1,nseg
if(segarea(jseg).le.zero) cycle
dij=dist(
I xyzseg(1,iseg),xyzseg(2,iseg),xyzseg(3,iseg),
J xyzseg(1,jseg),xyzseg(2,jseg),xyzseg(3,jseg))
if(dij.lt.dmin) then
c Accumulate the area at the centroid.
asum = segarea(iseg) + segarea(jseg)
do m=1,3
xyzseg(m,iseg) = (segarea(iseg)*xyzseg(m,iseg)+
+ segarea(jseg)*xyzseg(m,jseg))/asum
end do
segarea(iseg) = asum;
segarea(jseg) = zero;
c Start over, perhaps 3 or more segments should be merged
iseg = 0
exit
end if
end do
end do
endif
c Add the segments sorted by atom.
ief = 0
do iat=1,nat
do iseg=1,nseg
if(iat.ne.segatom(iseg).or.segarea(iseg).le.zero) cycle
c
ief = ief+1
dbl_mb(k_efcc+3*(ief-1)) =xyzseg(1,iseg)
dbl_mb(k_efcc+3*(ief-1)+1)=xyzseg(2,iseg)
dbl_mb(k_efcc+3*(ief-1)+2)=xyzseg(3,iseg)
dbl_mb(k_efcs+ief-1) = segarea(iseg) ! KS area
if (do_cosmo_model.eq.DO_COSMO_YK) then
c Store the YK variables.
zetai=zeta/sqrt(segarea(iseg)) ! zetai=zeta/sqrt(dsurf)/ratm_real
zetaii=zetai*sqrt(2d0/pi)
xyzff(iseg,iat) = one ! no segments inside other spheres
dbl_mb(k_efcs+ief-1) = zetaii ! zetaii/xyzff(iseg,iat))
dbl_mb(k_msrf+ief-1) = segarea(iseg) ! YK area
dbl_mb(k_efczz+ief-1) = zetai
end if
c Atom map.
int_mb(k_efciat+ief-1) = iat
enddo
enddo
nefc = ief
nseg = nefc
else
srfmol=0d0
volmol=0d0
ief =0
@ -805,6 +918,24 @@ c 775
enddo
endif
enddo
c Volume is underestimated with the above method, let's use GEPOL
c just to get a better volume estimate for the vdW cavity.
c For the KS method the SES cavity should be used.
if (do_cosmo_model.eq.DO_COSMO_YK) then
rd = 0! no solvent
rmim = rd/2! dummy, not used
ofac = 0.9 ! dummy, not used
dvec = 1.0 ! segment vector direction, not relevant
ndiv = 1
call gepol('WSURF',.false.,.false.,.false.,.false.,'TIT',0,
& rmim,'','','',ofac,rd,ndiv,dvec,
& .true.,.false.,nat,ratm,xyzatm,mxface,
& gseg,xyzseg,segarea,segatom,
& garea,gvolume,.false.)
volmol = gvolume
endif
c
endif
srfmol=srfmol*(cau2ang**2)
volmol=volmol*(cau2ang**3)
c
@ -928,7 +1059,8 @@ c
& errquit('cosmo_cossas chop stack k_efcc failed',911,MA_ERR)
c
return
9999 format(/,1x,'solvent accessible surface',/,1x,26(1h-))
9999 format(/,1x,' van der Waals surface ',/,1x,26(1h-))
8888 format(/,1x,' solvent excluded surface ',/,1x,26(1h-))
9998 format(/,1x,'---------- ATOMIC COORDINATES (A.U.) ----------',
1 '-- VDWR(ANG.) --')
9997 format( 1x,i5,3f14.8,f10.3)
@ -1794,13 +1926,13 @@ c
c
double precision allefc !< [Output] The total solute charge
!< density - COSMO charge interaction
!< energy
!< energy (pot*q)
double precision atmefc !< [Output] The nuclear - COSMO charge
!< interaction energy
!< interaction energy (Q*B*q)
double precision elcefc !< [Output] The electron - COSMO charge
!< interaction energy
double precision efcefc !< [Output] The COSMO charge - COSMO
!< charge interaction energy
!< charge interaction energy (q'*A*q/2)
c
c
integer ief, jef !< Counters over COSMO charges

View file

@ -40,20 +40,24 @@ c
integer ibegin,iend,ilen,i
character*4096 newchr
integer maxno_authors,len_newchr
parameter(maxno_authors=107)
parameter(maxno_authors=108)
character*64 name_author(maxno_authors)
data name_author/
+ 'E. Apra', 'E. J. Bylaska',
+ 'N. Govind', 'K. Kowalski',
+ 'D. Mejia-Rodriguez',
+ 'N. P. Bauman', 'A. Panyala',
+ 'R. J. Harrison',
+ 'M. Valiev',
+ 'D. Mejia-Rodriguez', 'A. Kunitsa',
+ 'N. P. Bauman', 'A. Panyala',
+ 'W. A. de Jong', 'T. P. Straatsma', 'H. J. J. van Dam',
+ 'D. Wang', 'T. L. Windus',
+ 'J. Hammond', 'J. Autschbach','A. Woods',
+ 'J. Hammond', 'J. Autschbach',
+ 'R. de P. Soares',
+ 'A. Kunitsa',
+ 'K. Bhaskaran-Nair',
+ 'J. Brabec','K. Lopata','S. A. Fischer','S. Krishnamoorthy',
+ 'M. Jacquelin',
+ 'A. Woods',
+ 'W. Ma', 'M. Klemm', 'O. Villa', 'Y. Chen', 'V. Anisimov',
+ 'F. Aquino', 'S. Hirata', 'M. T. Hackler',
+ 'E. Hermes', 'L. Jensen', 'J. E. Moore', 'J. C. Becca',
@ -61,7 +65,7 @@ c
+ 'M. Malagoli', 'A. Marenich', 'A. Otero-de-la-Roza',
+ 'J. Mullin', 'P. Nichols', 'R. Peverati', 'J. Pittner',
+ 'Y. Zhao', 'P.-D. Fan', 'A. Fonari', 'M. J. Williamson',
+ 'R. J. Harrison', 'J. R. Rehr',
+ 'J. R. Rehr',
+ 'M. Dupuis', 'D. Silverstein', 'D. M. A. Smith',
+ 'J. Nieplocha', 'V. Tipparaju', 'M. Krishnan',
+ 'B. E. Van Kuiken', 'A. Vazquez-Mayagoitia',

View file

@ -124,6 +124,7 @@ fi
fi
if [[ ! $(grep -i cosmo $TRAVIS_BUILD_DIR/src/stubs.F| awk '/cosmo_input/') ]]; then
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs cosmo_h2o_dft
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs cosmo_h2cco2_opt
fi
if [[ ! $(grep -i gw $TRAVIS_BUILD_DIR/src/stubs.F| awk '/gw_input/') ]]; then
cd $TRAVIS_BUILD_DIR/QA && ./runtests.mpi.unix procs $nprocs ritddft_h2o ritddft_co