followup of https://github.com/nwchemgit/nwchem/issues/776. Separation of SMD specific parameters

This commit is contained in:
edoapra 2023-08-09 18:11:59 -07:00
parent b8a3c97789
commit 1ee94b2d2b
No known key found for this signature in database
GPG key ID: 48E12DA1EDE9E1B0
3 changed files with 92 additions and 62 deletions

View file

@ -168,7 +168,7 @@ c --------------------------
> nmap*16*ma_sizeof(mt_byte,1,mt_byte),0)
call ga_sync()
if(.not.rtdb_put(rtdb,'cosmo:crad',mt_dbl, nmap,crad))
if(.not.rtdb_put(rtdb,'cosmo:crad',mt_dbl, nmap,crad(1)))
$ call errquit('cosmo_initialize: rtdb put crad ',914,
& rtdb_err)
c
@ -252,7 +252,8 @@ c
& errquit('cosmo_init,ma_pop_stack of l_tag failed',911, MA_ERR)
c
do_cosmo_model = DO_COSMO_YK
dielec =78.4d+00 ! default water dielectric
if(.not.rtdb_get(rtdb,'cosmo:dielec',mt_dbl,1,dielec))
D dielec =78.4d+00 ! default water dielectric
cosmo_screen = COSMO_SCREEN_ST
cosmo_sccor = COSMO_SCCOR_LAGRA ! surface charge correction adapted for geometry optimzations
minbem =2
@ -273,32 +274,32 @@ c MN solvation models -->
c
c reset some defaults if do_cosmo_smd is true
c
do_cosmo_smd = .false.
sola = -1.d0
solb = -1.d0
solc = -1.d0
solg = -1.d0
solh = -1.d0
soln = -1.d0
if (dielec.eq.78.4d+00) then
solvname_short = 'h2o'
solvname_long = 'water'
sola = 0.82d0
solb = 0.35d0
solc = 0.d0
solg = 103.62d0
solh = 0.d0
soln = 1.333d0
dielecinf = soln*soln ! n**2 (refractive index squarred)
end if
status =
$ rtdb_get(rtdb,'cosmo:do_cosmo_smd',mt_log,1,do_cosmo_smd)
c
if(.not. rtdb_get(rtdb,
O 'cosmo:do_cosmo_smd',mt_log,1,do_cosmo_smd))
D do_cosmo_smd=.false.
if (do_cosmo_smd) then
minbem =3
maxbem =4
ificos =0
end if
minbem =3
maxbem =4
ificos =0
if (dielec.eq.78.4d+00) then
solvname_short = 'h2o'
solvname_long = 'water'
sola = 0.82d0
solb = 0.35d0
solc = 0.d0
solg = 103.62d0
solh = 0.d0
soln = 1.333d0
dielecinf = soln*soln ! n**2 (refractive index squared)
else
sola = -1.d0
solb = -1.d0
solc = -1.d0
solg = -1.d0
solh = -1.d0
soln = -1.d0
end if
endif
c
status = rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo)
status = rtdb_get(rtdb,'cosmo:rsolv',mt_dbl,1,rsolv)
@ -320,7 +321,10 @@ c
status = rtdb_get(rtdb,'cosmo:zeta',mt_dbl,1,zeta)
status = rtdb_get(rtdb,'cosmo:gamma_s',mt_dbl,1,gammas)
status = rtdb_get(rtdb,'cosmo:sw_tol',mt_dbl,1,swtol)
status = rtdb_cget(rtdb,'cosmo:solvent',1,solvname_short)
if(.not.rtdb_cget(rtdb,'cosmo:solvent',1,solvname_short))
c solvname_short='unkn'
if(.not.rtdb_cget(rtdb,'cosmo:solvname_long',1,solvname_long))
$ solvname_long='unknown'
status = rtdb_get(rtdb,'cosmo:iangleb',mt_int,1,iangleb)
status = rtdb_get(rtdb,'cosmo:fibonacci',mt_int,1,fibonacci)
status = rtdb_get(rtdb,'cosmo:thomson',mt_int,1,thomson)
@ -341,37 +345,50 @@ c
c if set by the user
c
status = rtdb_get(rtdb,'cosmo:dielec',mt_dbl,1,dielec)
status = status.and.
S rtdb_get(rtdb,'cosmo:dielecinf',mt_dbl,1,dielecinf)
status = status.and.
S rtdb_get(rtdb,'cosmo:sola',mt_dbl,1,sola)
status = status.and.
S rtdb_get(rtdb,'cosmo:solb',mt_dbl,1,solb)
status = status.and.
S rtdb_get(rtdb,'cosmo:solc',mt_dbl,1,solc)
status = status.and.
S rtdb_get(rtdb,'cosmo:solg',mt_dbl,1,solg)
status = status.and.
S rtdb_get(rtdb,'cosmo:solh',mt_dbl,1,solh)
status = status.and.
S rtdb_get(rtdb,'cosmo:soln',mt_dbl,1,soln)
status = status.and.
S rtdb_get(rtdb,'cosmo:icds',mt_int,1,icds)
if(.not.status) then
if(odbug.and.me.eq.0)
& write(luout,*) "calling solv_data"
if(.not.solv_data(solvname_short,solvname_long, !short and long solvent names
& dielec,dielecinf, ! slow and fast components
& sola,solb,solc,solg,solh,soln,icds)) then
call errquit('solv_data failure',0,0)
if(do_cosmo_smd) then
status = status.or.
S rtdb_get(rtdb,'cosmo:dielecinf',mt_dbl,1,dielecinf)
status = status.or.
S rtdb_get(rtdb,'cosmo:sola',mt_dbl,1,sola)
status = status.or.
S rtdb_get(rtdb,'cosmo:solb',mt_dbl,1,solb)
status = status.or.
S rtdb_get(rtdb,'cosmo:solc',mt_dbl,1,solc)
status = status.or.
S rtdb_get(rtdb,'cosmo:solg',mt_dbl,1,solg)
status = status.or.
S rtdb_get(rtdb,'cosmo:solh',mt_dbl,1,solh)
status = status.or.
S rtdb_get(rtdb,'cosmo:soln',mt_dbl,1,soln)
#if 1
if(rtdb_get(rtdb,'cosmo:icds',mt_int,1,icds)) then
status=.true.
endif
#endif
endif
c status = status.or.
c S rtdb_get(rtdb,'cosmo:icds',mt_int,1,icds)
if (solvname_short.eq.'unkn') then
solvname_long='unknown'
else
if(.not.status) then
if(odbug.and.me.eq.0)
& write(luout,*) "calling solv_data for solv:",
S solvname_short,'@',solvname_long,'@'
if(.not.solv_data(solvname_short,solvname_long, !short and long solvent names
& dielec,dielecinf, ! slow and fast components
& sola,solb,solc,solg,solh,soln,icds)) then
call errquit('solv_data failure',0,0)
endif
endif
endif
c
c set sola,solb,solc,solg,solh,soln,icds parameters (either from
c solv_data or user-defined)
c
if (.not. rtdb_put(rtdb,'cosmo:dielec',mt_dbl,1,dielec))
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
if(do_cosmo_smd) then
if (.not. rtdb_put(rtdb,'cosmo:dielecinf',mt_dbl,1,dielecinf))
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
if (.not. rtdb_put(rtdb,'cosmo:sola',mt_dbl,1,sola))
@ -387,7 +404,10 @@ c
if (.not. rtdb_put(rtdb,'cosmo:soln',mt_dbl,1,soln))
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
if (.not. rtdb_put(rtdb,'cosmo:icds',mt_int,1,icds))
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
$ call errquit('cosmo_initialize: rtdb_put failed',0, RTDB_ERR)
if (.not.rtdb_cput(rtdb,'cosmo:solvname_long',1,solvname_long))
$ call errquit('cosmo_initialize:rtdb_cput failed',0, RTDB_ERR)
endif
c
if (me.eq.0) then
c write(luout,'(A)') " after solv_data"
@ -401,8 +421,8 @@ c write(luout,'(A)') " after solv_data"
if (dielecinf.lt.zero) call errquit
$ ('cosmo_initialize: illegal/undefined value for dielecinf',0,0)
write(luout,'(A,F8.4)') " dielec: ",dielec
write(luout,'(A,F8.4)') " dielecinf: ",dielecinf
if (do_cosmo_smd) then
write(luout,'(A,F8.4)') " dielecinf: ",dielecinf
if (icds.eq.1) then
write(luout,'(A)') " aqueous SMD model solvent descriptors"
write(luout,'(A,F8.4)') " dielec: ",dielec
@ -433,7 +453,7 @@ c write(luout,'(A)') " after solv_data"
T then
icds = 2
else
write(luout,'(A,I1)') " icds: ",icds
write(luout,'(a,i1)') " icds: ",icds
call errquit
$ ('cosmo_initialize: illegal value for icds',0,0)
endif

View file

@ -593,13 +593,19 @@ c
c
oprint_energies = util_print("cosmo_energies",print_never)
oprint_energies = oprint_energies.and.(ga_nodeid().eq.0)
c
wgsrf_cosmo_vem = zero
wstar_cosmo_vem = zero
tolw = 0.00001d0
iter_cosmo_vem = 0
dvem1=(dielecinf-one)*dielec/dielecinf/(dielec-one)
dvem2=(dielec-dielecinf)/dielecinf/(dielec-one)
if(dielecinf.ne.0d0) then
dvem1=(dielecinf-1d0)*dielec/dielecinf/(dielec-1d0)
dvem2=(dielec-dielecinf)/dielecinf/(dielec-1d0)
else
dvem1=0d0
dvem2=0d0
endif
c
if (.not. rtdb_get
$ (rtdb,'cosmo:istep_cosmo_vem',mt_int,1,istep_cosmo_vem))
@ -1263,8 +1269,10 @@ c calculated in tddft_grad_compute_g.F.
c
c calculate noneq charges from current ES charges and saved GS charges
c
if (do_cosmo_vem.eq.2) dvem1=one
if (do_cosmo_vem.eq.2) dvem2=zero
if (do_cosmo_vem.eq.2) then
dvem1=1d0
dvem2=0d0
endif
espolgsrf = zero
espolin = zero
espoldyn = zero

View file

@ -17,8 +17,12 @@ c
c
c get the solvent data
c
write(6,*)' SSSSSSSSSSSSSSS'
if (solvname_short.eq.'h2o') solvname_long='water' ! default
if (solvname_short.eq.'unkn') then
solvname_long='unknown' ! default
solv_data=.false.
return
endif
eps=0d0
epsinf=0d0
call DATSOLMN(solvname_short,solvname_long,
@ -36,8 +40,6 @@ c
if (solvname_short.ne.'h2o'.and.solvname_long.ne.'water')
& icds = 2
c
write(6,*) ' icds ',icds
write(6,*) solvname_short, solvname_long
return
end
c