mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-27 13:45:27 -04:00
added c0 parameter to penalty function. Was Vpenalty = K0*exp(-((gamm-gamm0)/gcut0)**2) now is K0*exp(-((gamm-gamm0)/gcut0)**c0)...EJB
This commit is contained in:
Eric Bylaska
parent
3ce377d4fc
commit
65810a6726
8 changed files with 171 additions and 51 deletions
|
|
@ -768,9 +768,9 @@ c ---------------
|
|||
integer inb,nhb,nhc,nhp
|
||||
integer iat,jat
|
||||
double precision r
|
||||
double precision r0,k,K0,gcut0,gamma0,gamm
|
||||
double precision r0,k,K0,gcut0,gamma0,c0,gamm
|
||||
integer i_c,i_rhb,i_khb,i_ijhb
|
||||
integer i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0
|
||||
integer i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0,i_c0
|
||||
double precision etot,e,f(3)
|
||||
character*255 namespace
|
||||
logical cons_active
|
||||
|
|
@ -883,6 +883,7 @@ c
|
|||
i_k0 = cons_get_i_pbondings_k0()
|
||||
i_gcut0 = cons_get_i_pbondings_gcut0()
|
||||
i_gamma0 = cons_get_i_pbondings_gamma0()
|
||||
i_c0 = cons_get_i_pbondings_c0()
|
||||
|
||||
ii = 0
|
||||
do i=1,nhp
|
||||
|
|
@ -890,14 +891,16 @@ c
|
|||
K0 = dbl_mb(i_k0+i-1)
|
||||
gcut0 = dbl_mb(i_gcut0+i-1)
|
||||
gamma0 = dbl_mb(i_gamma0+i-1)
|
||||
c0 = dbl_mb(i_c0+i-1)
|
||||
|
||||
call cons_pbondings_energy(n,int_mb(i_indx+2*ii),
|
||||
> dbl_mb(i_coef+ii),
|
||||
> K0,gcut0,gamma0,dbl_mb(i_c),gamm,e)
|
||||
> K0,gcut0,gamma0,c0,
|
||||
> dbl_mb(i_c),gamm,e)
|
||||
if(master) then
|
||||
write(luout,'(A,I5)')" bondings penalty# :",i
|
||||
write(luout,'(A,3(F12.6,2X))')
|
||||
> " spring parameters (K0,gcut0,gamma0):",K0,gcut0,gamma0
|
||||
write(luout,'(A,4(F12.6,2X))')
|
||||
> " spring parameters (K0,gcut0,gamma0,c0):",K0,gcut0,gamma0,c0
|
||||
write(luout,'(A,F12.6)')
|
||||
> " gamma :",gamm
|
||||
write(luout,'(A,F12.6)')
|
||||
|
|
|
|||
|
|
@ -98,6 +98,11 @@
|
|||
integer cons_get_h_pbondings_gamma0
|
||||
external cons_get_h_pbondings_gamma0
|
||||
|
||||
integer cons_get_i_pbondings_c0
|
||||
external cons_get_i_pbondings_c0
|
||||
integer cons_get_h_pbondings_c0
|
||||
external cons_get_h_pbondings_c0
|
||||
|
||||
|
||||
|
||||
logical cons_fixed_atoms
|
||||
|
|
|
|||
|
|
@ -623,8 +623,6 @@ c cons_get_i_c = i_c
|
|||
end
|
||||
|
||||
|
||||
|
||||
|
||||
integer function cons_get_h_pbondings_gamma0()
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
|
|
@ -644,6 +642,29 @@ c cons_get_i_c = i_c
|
|||
return
|
||||
end
|
||||
|
||||
|
||||
integer function cons_get_h_pbondings_c0()
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
cons_get_h_pbondings_c0 = h_pbondings_c0
|
||||
return
|
||||
end
|
||||
|
||||
integer function cons_get_i_pbondings_c0()
|
||||
implicit none
|
||||
#include "mafdecls.fh"
|
||||
#include "cons_data.fh"
|
||||
#include "errquit.fh"
|
||||
if(.not. ma_get_index(h_pbondings_c0,
|
||||
> cons_get_i_pbondings_c0))
|
||||
> call errquit('cons_data: unable to get i_pbondings_c0',
|
||||
> h_pbondings_c0, RTDB_ERR)
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine cons_set_pbondings_nhp(n)
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
|
|
@ -742,7 +763,6 @@ c cons_get_i_c = i_c
|
|||
end
|
||||
|
||||
|
||||
|
||||
subroutine cons_set_pbondings_gamma0(n)
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
|
|
@ -760,6 +780,22 @@ c cons_get_i_c = i_c
|
|||
end
|
||||
|
||||
|
||||
subroutine cons_set_pbondings_c0(n)
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
integer n
|
||||
h_pbondings_c0=n
|
||||
return
|
||||
end
|
||||
|
||||
subroutine cons_get_pbondings_c0(n)
|
||||
implicit none
|
||||
#include "cons_data.fh"
|
||||
integer n
|
||||
n=h_pbondings_c0
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -43,6 +43,7 @@ c parameter (max_fixed_atoms=100)
|
|||
integer h_pbondings_k0
|
||||
integer h_pbondings_gcut0
|
||||
integer h_pbondings_gamma0
|
||||
integer h_pbondings_c0
|
||||
|
||||
real*8 unita(3,3),ub(3,3)
|
||||
integer isystype
|
||||
|
|
@ -61,7 +62,7 @@ c parameter (max_fixed_atoms=100)
|
|||
& h_pbondings_n0,h_pbondings_indx,
|
||||
& h_pbondings_coef,
|
||||
& h_pbondings_k0,h_pbondings_gcut0,
|
||||
& h_pbondings_gamma0,
|
||||
& h_pbondings_gamma0,h_pbondings_c0,
|
||||
& h_fpos,i_fpos,
|
||||
& h_c,
|
||||
& cons_irtdb,
|
||||
|
|
|
|||
|
|
@ -19,7 +19,7 @@ c
|
|||
integer nfp
|
||||
logical status
|
||||
integer id(4)
|
||||
double precision k,r,K0,gcut0,gamma0,coef(20)
|
||||
double precision k,r,K0,gcut0,gamma0,c0,coef(20)
|
||||
integer n,indx(20)
|
||||
integer geom, ncenters
|
||||
c
|
||||
|
|
@ -176,10 +176,11 @@ c ----------------------------
|
|||
> call errquit('cons_input: failed to read keyword', 0,
|
||||
> 0)
|
||||
if (inp_compare(.false.,'bondings', test)) then
|
||||
if(cons_penalty_bondings_input(n,indx,coef,K0,gcut0,gamma0))
|
||||
if(cons_penalty_bondings_input(n,indx,coef,
|
||||
> K0,gcut0,gamma0,c0))
|
||||
> then
|
||||
call cons_add_pbondings(namespace,n,indx,coef,K0,
|
||||
> gcut0,gamma0,rtdb)
|
||||
> gcut0,gamma0,c0,rtdb)
|
||||
else
|
||||
call errquit('cons_input:failed penalty bondings input',0,0)
|
||||
end if
|
||||
|
|
@ -371,10 +372,10 @@ c
|
|||
* * *
|
||||
* ****************************************
|
||||
*
|
||||
* pbondings K0 gcut0 gamma0 -1.0 1 2 1.0 2 3
|
||||
* penalty bondings K0 gcut0 gamma0 c0 -1.0 1 2 1.0 2 3
|
||||
*
|
||||
logical function cons_penalty_bondings_input(n,indx,coef,K0,
|
||||
> gcut0,gamma0)
|
||||
> gcut0,gamma0,c0)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "inp.fh"
|
||||
|
|
@ -384,20 +385,22 @@ c
|
|||
c
|
||||
logical ok
|
||||
integer n,indx(*)
|
||||
double precision coef(*),K0,gcut0,gamma0
|
||||
double precision coef(*),K0,gcut0,gamma0,c0
|
||||
|
||||
integer itest
|
||||
double precision ftest
|
||||
|
||||
ok = .true.
|
||||
|
||||
!*** enter K0, gcut0, and gamma0 ***
|
||||
!*** enter K0, gcut0, gamma0, and c0 ***
|
||||
if (.not.inp_f(ftest)) ok = .false.
|
||||
K0 = ftest
|
||||
if (ok.and.(.not.inp_f(ftest))) ok = .false.
|
||||
gcut0 = ftest
|
||||
if (ok.and.(.not.inp_f(ftest))) ok = .false.
|
||||
gamma0 = ftest
|
||||
if (ok.and.(.not.inp_f(ftest))) ok = .false.
|
||||
c0 = ftest
|
||||
|
||||
!*** coef and indx ***
|
||||
n = 0
|
||||
|
|
|
|||
|
|
@ -1312,7 +1312,8 @@ c
|
|||
integer h_tmp_k0,i_tmp_k0
|
||||
integer h_tmp_gcut0,i_tmp_gcut0
|
||||
integer h_tmp_gamma0,i_tmp_gamma0
|
||||
double precision K0, gcut0,gamma0
|
||||
integer h_tmp_c0,i_tmp_c0
|
||||
double precision K0, gcut0,gamma0,c0
|
||||
character*255 tag_n
|
||||
character*255 tag_n0
|
||||
character*255 tag_indx
|
||||
|
|
@ -1320,6 +1321,7 @@ c
|
|||
character*255 tag_k0
|
||||
character*255 tag_gcut0
|
||||
character*255 tag_gamma0
|
||||
character*255 tag_c0
|
||||
|
||||
call cons_pbondings_n_tag(namespace,tag_n)
|
||||
call cons_pbondings_n0_tag(namespace,tag_n0)
|
||||
|
|
@ -1328,6 +1330,7 @@ c
|
|||
call cons_pbondings_k0_tag(namespace,tag_k0)
|
||||
call cons_pbondings_gcut0_tag(namespace,tag_gcut0)
|
||||
call cons_pbondings_gamma0_tag(namespace,tag_gamma0)
|
||||
call cons_pbondings_c0_tag(namespace,tag_c0)
|
||||
|
||||
|
||||
if(.not.rtdb_get(rtdb,tag_n,mt_int,1,nhp)) then
|
||||
|
|
@ -1374,6 +1377,11 @@ c
|
|||
if ( .not. ma_push_get( MT_DBL, nhp, 'h_tmp_gamma0',
|
||||
& h_tmp_gamma0, i_tmp_gamma0) ) call errquit(
|
||||
& 'cons_info_pbondings: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
|
||||
if ( .not. ma_push_get( MT_DBL, nhp, 'h_tmp_c0',
|
||||
& h_tmp_c0, i_tmp_c0) ) call errquit(
|
||||
& 'cons_info_pbondings: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
|
||||
status = rtdb_get(rtdb,tag_indx,mt_int,2*mm,int_mb(i_tmp_indx))
|
||||
|
|
@ -1385,6 +1393,9 @@ c
|
|||
> rtdb_get(rtdb,tag_gcut0,mt_dbl,nhp,dbl_mb(i_tmp_gcut0))
|
||||
status = status.and.
|
||||
> rtdb_get(rtdb,tag_gamma0,mt_dbl,nhp,dbl_mb(i_tmp_gamma0))
|
||||
status = status.and.
|
||||
> rtdb_get(rtdb,tag_c0,mt_dbl,nhp,dbl_mb(i_tmp_c0))
|
||||
|
||||
if(.not.status) call errquit(
|
||||
> 'cons_info_pbondings: unable to get prior harm bondings',
|
||||
> nhp, MA_ERR)
|
||||
|
|
@ -1396,12 +1407,13 @@ c
|
|||
K0 = dbl_mb(i_tmp_k0+i-1)
|
||||
gcut0 = dbl_mb(i_tmp_gcut0+i-1)
|
||||
gamma0 = dbl_mb(i_tmp_gamma0+i-1)
|
||||
c0 = dbl_mb(i_tmp_c0+i-1)
|
||||
|
||||
write(luout,*)
|
||||
write(luout,'(A,I5)') " Bondings Penalty #",i
|
||||
write(luout,'(A,F18.9,F18.9,E14.6)')
|
||||
> " - parameters (K0,gcut0,gamma0):",
|
||||
> K0,gcut0,gamma0
|
||||
write(luout,'(A,F18.9,F18.9,E14.6,F18.9)')
|
||||
> " - parameters (K0,gcut0,gamma0,c0):",
|
||||
> K0,gcut0,gamma0,c0
|
||||
write(luout,*) " - coefficient index1 index2"
|
||||
do j=1,n
|
||||
write(luout,'(F16.6,I7,I7)') dbl_mb(i_tmp_coef+ii+j-1),
|
||||
|
|
@ -1413,6 +1425,9 @@ c
|
|||
end do
|
||||
end if
|
||||
|
||||
if (.not. ma_pop_stack(h_tmp_c0) ) call errquit(
|
||||
& 'cons_info_pbondings: unable to free h_tmp_c0',
|
||||
& 0, MA_ERR)
|
||||
if (.not. ma_pop_stack(h_tmp_gamma0) ) call errquit(
|
||||
& 'cons_info_pbondings: unable to free h_tmp_gamma0',
|
||||
& 0, MA_ERR)
|
||||
|
|
@ -1444,11 +1459,11 @@ c
|
|||
* ************************************
|
||||
|
||||
subroutine cons_add_pbondings(namespace,n,indx,coef,
|
||||
> K0,gcut0,gamma0,rtdb)
|
||||
> K0,gcut0,gamma0,c0,rtdb)
|
||||
implicit none
|
||||
character*(*) namespace
|
||||
integer n,indx(*)
|
||||
double precision coef(*),K0,gcut0,gamma0
|
||||
double precision coef(*),K0,gcut0,gamma0,c0
|
||||
integer rtdb
|
||||
|
||||
#include "errquit.fh"
|
||||
|
|
@ -1467,6 +1482,7 @@ c
|
|||
integer h_tmp_k0,i_tmp_k0
|
||||
integer h_tmp_gcut0,i_tmp_gcut0
|
||||
integer h_tmp_gamma0,i_tmp_gamma0
|
||||
integer h_tmp_c0,i_tmp_c0
|
||||
character*255 tag_n
|
||||
character*255 tag_n0
|
||||
character*255 tag_indx
|
||||
|
|
@ -1474,6 +1490,7 @@ c
|
|||
character*255 tag_k0
|
||||
character*255 tag_gcut0
|
||||
character*255 tag_gamma0
|
||||
character*255 tag_c0
|
||||
|
||||
call cons_pbondings_n_tag(namespace,tag_n)
|
||||
call cons_pbondings_n0_tag(namespace,tag_n0)
|
||||
|
|
@ -1482,6 +1499,7 @@ c
|
|||
call cons_pbondings_k0_tag(namespace,tag_k0)
|
||||
call cons_pbondings_gcut0_tag(namespace,tag_gcut0)
|
||||
call cons_pbondings_gamma0_tag(namespace,tag_gamma0)
|
||||
call cons_pbondings_c0_tag(namespace,tag_c0)
|
||||
|
||||
|
||||
if(.not.rtdb_get(rtdb,tag_n,mt_int,1,nhp0)) then
|
||||
|
|
@ -1514,6 +1532,10 @@ c
|
|||
if (.not. ma_push_get( MT_DBL, nhp, 'h_tmp_gamma0',
|
||||
& h_tmp_gamma0, i_tmp_gamma0) ) call errquit(
|
||||
& 'cons_add_spring: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
if (.not. ma_push_get( MT_DBL, nhp, 'h_tmp_c0',
|
||||
& h_tmp_c0, i_tmp_c0) ) call errquit(
|
||||
& 'cons_add_spring: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
|
||||
else
|
||||
|
|
@ -1545,6 +1567,10 @@ c
|
|||
if ( .not. ma_push_get( MT_DBL, nhp, 'h_tmp_gamma0',
|
||||
& h_tmp_gamma0, i_tmp_gamma0) ) call errquit(
|
||||
& 'cons_add_spring: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
if ( .not. ma_push_get( MT_DBL, nhp, 'h_tmp_c0',
|
||||
& h_tmp_c0, i_tmp_c0) ) call errquit(
|
||||
& 'cons_add_spring: unable to allocate scratch space',
|
||||
& nhp, MA_ERR)
|
||||
|
||||
status = status . and. rtdb_get(rtdb,tag_indx,
|
||||
|
|
@ -1557,6 +1583,9 @@ c
|
|||
> mt_dbl,nhp0,dbl_mb(i_tmp_gcut0))
|
||||
status = status . and. rtdb_get(rtdb,tag_gamma0,
|
||||
> mt_dbl,nhp0,dbl_mb(i_tmp_gamma0))
|
||||
status = status . and. rtdb_get(rtdb,tag_c0,
|
||||
> mt_dbl,nhp0,dbl_mb(i_tmp_c0))
|
||||
|
||||
if(.not.status) call errquit(
|
||||
> 'cons_add_spring: unable to get prior harm bondings',
|
||||
> nhp, MA_ERR)
|
||||
|
|
@ -1568,6 +1597,7 @@ c
|
|||
status = status .and. rtdb_delete(rtdb,tag_k0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_gcut0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_gamma0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_c0)
|
||||
if(.not.status) call errquit(
|
||||
> 'cons_add_spring: unable to delete prior harm bondingss',
|
||||
> 0, MA_ERR)
|
||||
|
|
@ -1597,11 +1627,16 @@ c
|
|||
> mt_dbl,nhp,dbl_mb(i_tmp_gcut0))
|
||||
status = status.and. rtdb_put(rtdb,tag_gamma0,
|
||||
> mt_dbl,nhp,dbl_mb(i_tmp_gamma0))
|
||||
status = status.and. rtdb_put(rtdb,tag_c0,
|
||||
> mt_dbl,nhp,dbl_mb(i_tmp_c0))
|
||||
|
||||
if(.not.status) call errquit(
|
||||
> 'cons_add_spring: unable to get prior harm bondings',
|
||||
> nhp, MA_ERR)
|
||||
|
||||
if (.not. ma_pop_stack(h_tmp_c0) ) call errquit(
|
||||
& 'cons_add_spring: unable to free h_tmp_c0',
|
||||
& 0, MA_ERR)
|
||||
if (.not. ma_pop_stack(h_tmp_gamma0) ) call errquit(
|
||||
& 'cons_add_spring: unable to free h_tmp_gamma0',
|
||||
& 0, MA_ERR)
|
||||
|
|
|
|||
|
|
@ -542,6 +542,7 @@ c
|
|||
integer h_pbondings_k0
|
||||
integer h_pbondings_gcut0
|
||||
integer h_pbondings_gamma0
|
||||
integer h_pbondings_c0
|
||||
|
||||
h_pbondings_n0 = cons_get_h_pbondings_n0()
|
||||
h_pbondings_indx = cons_get_h_pbondings_indx()
|
||||
|
|
@ -549,6 +550,7 @@ c
|
|||
h_pbondings_k0 = cons_get_h_pbondings_k0()
|
||||
h_pbondings_gcut0 = cons_get_h_pbondings_gcut0()
|
||||
h_pbondings_gamma0 = cons_get_h_pbondings_gamma0()
|
||||
h_pbondings_c0 = cons_get_h_pbondings_c0()
|
||||
c
|
||||
c unload harmonic constraints if any
|
||||
c
|
||||
|
|
@ -571,6 +573,9 @@ c
|
|||
& 0, MA_ERR)
|
||||
if (.not. ma_free_heap(h_pbondings_gamma0) ) call errquit(
|
||||
& 'cons_bonding_input: unable to free h_pbondings_gamma0',
|
||||
& 0, MA_ERR)
|
||||
if (.not. ma_free_heap(h_pbondings_c0) ) call errquit(
|
||||
& 'cons_bonding_input: unable to free h_pbondings_c0',
|
||||
& 0, MA_ERR)
|
||||
call cons_set_pbondings_nhp(0)
|
||||
end if
|
||||
|
|
@ -596,6 +601,7 @@ c
|
|||
integer h_pbondings_k0,i_pbondings_k0
|
||||
integer h_pbondings_gcut0,i_pbondings_gcut0
|
||||
integer h_pbondings_gamma0,i_pbondings_gamma0
|
||||
integer h_pbondings_c0,i_pbondings_c0
|
||||
character*255 tag_n
|
||||
character*255 tag_n0
|
||||
character*255 tag_indx
|
||||
|
|
@ -603,6 +609,7 @@ c
|
|||
character*255 tag_k0
|
||||
character*255 tag_gcut0
|
||||
character*255 tag_gamma0
|
||||
character*255 tag_c0
|
||||
integer ma_type,ma_n
|
||||
c
|
||||
call cons_unload_pbondings()
|
||||
|
|
@ -615,6 +622,7 @@ c
|
|||
call cons_pbondings_k0_tag(namespace,tag_k0)
|
||||
call cons_pbondings_gcut0_tag(namespace,tag_gcut0)
|
||||
call cons_pbondings_gamma0_tag(namespace,tag_gamma0)
|
||||
call cons_pbondings_c0_tag(namespace,tag_c0)
|
||||
c
|
||||
c load harmonic constraints
|
||||
c
|
||||
|
|
@ -682,6 +690,14 @@ c > nhb, MA_ERR)
|
|||
& 'cons_load_hbonds: unable to allocate cons scratch space',
|
||||
& 2*nhp, MA_ERR)
|
||||
|
||||
if ( .not. rtdb_ma_get( rtdb, tag_c0,ma_type, ma_n,
|
||||
& h_pbondings_c0) ) call errquit(
|
||||
& 'cons_load_hbonds: unable to allocate cons scratch space',
|
||||
& 2*nhp, MA_ERR)
|
||||
if ( .not. ma_get_index(h_pbondings_c0,
|
||||
& i_pbondings_c0) ) call errquit(
|
||||
& 'cons_load_hbonds: unable to allocate cons scratch space',
|
||||
& 2*nhp, MA_ERR)
|
||||
|
||||
call cons_set_pbondings_nhp(nhp)
|
||||
call cons_set_pbondings_n0(h_pbondings_n0)
|
||||
|
|
@ -690,6 +706,7 @@ c > nhb, MA_ERR)
|
|||
call cons_set_pbondings_k0(h_pbondings_k0)
|
||||
call cons_set_pbondings_gcut0(h_pbondings_gcut0)
|
||||
call cons_set_pbondings_gamma0(h_pbondings_gamma0)
|
||||
call cons_set_pbondings_c0(h_pbondings_c0)
|
||||
|
||||
return
|
||||
end
|
||||
|
|
@ -711,6 +728,7 @@ c
|
|||
character*255 tag_k0
|
||||
character*255 tag_gcut0
|
||||
character*255 tag_gamma0
|
||||
character*255 tag_c0
|
||||
|
||||
call cons_pbondings_n_tag(namespace,tag_n)
|
||||
call cons_pbondings_n0_tag(namespace,tag_n0)
|
||||
|
|
@ -719,6 +737,7 @@ c
|
|||
call cons_pbondings_k0_tag(namespace,tag_k0)
|
||||
call cons_pbondings_gcut0_tag(namespace,tag_gcut0)
|
||||
call cons_pbondings_gamma0_tag(namespace,tag_gamma0)
|
||||
call cons_pbondings_c0_tag(namespace,tag_c0)
|
||||
|
||||
status = rtdb_delete(rtdb,tag_n)
|
||||
status = status .and. rtdb_delete(rtdb,tag_n0)
|
||||
|
|
@ -727,6 +746,7 @@ c
|
|||
status = status .and. rtdb_delete(rtdb,tag_k0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_gcut0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_gamma0)
|
||||
status = status .and. rtdb_delete(rtdb,tag_c0)
|
||||
|
||||
end
|
||||
|
||||
|
|
@ -746,9 +766,9 @@ c
|
|||
double precision energy
|
||||
logical oprint
|
||||
integer i,ii,n,j
|
||||
double precision K0,gcut0,gamma0,gamm
|
||||
double precision K0,gcut0,gamma0,c0,gamm
|
||||
double precision e,enb
|
||||
integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0
|
||||
integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0,i_c0
|
||||
integer nhp
|
||||
|
||||
c
|
||||
|
|
@ -765,6 +785,7 @@ c
|
|||
i_k0 = cons_get_i_pbondings_k0()
|
||||
i_gcut0 = cons_get_i_pbondings_gcut0()
|
||||
i_gamma0 = cons_get_i_pbondings_gamma0()
|
||||
i_c0 = cons_get_i_pbondings_c0()
|
||||
|
||||
e=0.0d0
|
||||
ii = 0
|
||||
|
|
@ -773,16 +794,17 @@ c
|
|||
K0 = dbl_mb(i_k0+i-1)
|
||||
gcut0 = dbl_mb(i_gcut0+i-1)
|
||||
gamma0 = dbl_mb(i_gamma0+i-1)
|
||||
c0 = dbl_mb(i_c0+i-1)
|
||||
|
||||
call cons_pbondings_energy(n,int_mb(i_indx+2*ii),
|
||||
> dbl_mb(i_coef+ii),
|
||||
> K0,gcut0,gamma0,dbl_mb(i_c),gamm,enb)
|
||||
> dbl_mb(i_coef+ii),
|
||||
> K0,gcut0,gamma0,c0,dbl_mb(i_c),gamm,enb)
|
||||
if(oprint) then
|
||||
write(luout,*)
|
||||
write(luout,'(A,I4)') "Adding bondings spring # ",i
|
||||
write(luout,'(A,F18.9,F18.9,E14.6)')
|
||||
> " - spring parameters (K0,gcut0,gamma0):",
|
||||
> K0,gcut0,gamma0
|
||||
write(luout,'(A,F18.9,F18.9,E14.6,F18.9)')
|
||||
> " - spring parameters (K0,gcut0,gamma0,c0):",
|
||||
> K0,gcut0,gamma0,c0
|
||||
write(luout,'(A,E18.6)')
|
||||
> " - gamma :",gamm
|
||||
write(luout,'(A,F18.9)')
|
||||
|
|
@ -828,9 +850,9 @@ c ---------------
|
|||
logical oprint
|
||||
integer i,j,n,ii
|
||||
integer nhp
|
||||
double precision gamma0,K0,gcut0,gamm
|
||||
double precision gamma0,K0,gcut0,c0,gamm
|
||||
double precision e,enb
|
||||
integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0
|
||||
integer i_c,i_n,i_indx,i_coef,i_k0,i_gcut0,i_gamma0,i_c0
|
||||
c
|
||||
call util_print_push()
|
||||
call util_print_rtdb_load(rtdb,'cons')
|
||||
|
|
@ -845,6 +867,7 @@ c
|
|||
i_k0 = cons_get_i_pbondings_k0()
|
||||
i_gcut0 = cons_get_i_pbondings_gcut0()
|
||||
i_gamma0 = cons_get_i_pbondings_gamma0()
|
||||
i_c0 = cons_get_i_pbondings_c0()
|
||||
|
||||
|
||||
e=0.0d0
|
||||
|
|
@ -854,16 +877,17 @@ c
|
|||
K0 = dbl_mb(i_k0+i-1)
|
||||
gcut0 = dbl_mb(i_gcut0+i-1)
|
||||
gamma0 = dbl_mb(i_gamma0+i-1)
|
||||
c0 = dbl_mb(i_c0+i-1)
|
||||
|
||||
call cons_pbondings_force(n,int_mb(i_indx+2*ii),
|
||||
> dbl_mb(i_coef+ii),
|
||||
> K0,gcut0,gamma0,dbl_mb(i_c),gamm,enb,gx)
|
||||
> dbl_mb(i_coef+ii),
|
||||
> K0,gcut0,gamma0,c0,dbl_mb(i_c),gamm,enb,gx)
|
||||
if(oprint) then
|
||||
write(luout,*)
|
||||
write(luout,'(A,I4)') "Adding bondings penalty # ",i
|
||||
write(luout,'(A,F18.9,F18.9,E14.6)')
|
||||
> " - spring parameters (K0,gcut0,gamma0):",
|
||||
> K0,gcut0,gamma0
|
||||
write(luout,'(A,F18.9,F18.9,E14.6,F18.9)')
|
||||
> " - spring parameters (K0,gcut0,gamma0,c0):",
|
||||
> K0,gcut0,gamma0,c0
|
||||
write(luout,'(A,E18.6)')
|
||||
> " - gamma :",gamm
|
||||
write(luout,'(A,F18.9)')
|
||||
|
|
@ -1287,8 +1311,8 @@ c
|
|||
|
||||
|
||||
|
||||
subroutine cons_pbondings_energy(n,indx,coef,K0,gcut0,gamma0,rion,
|
||||
> gamm,energy)
|
||||
subroutine cons_pbondings_energy(n,indx,coef,K0,gcut0,gamma0,c0,
|
||||
> rion,gamm,energy)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
|
|
@ -1297,7 +1321,7 @@ c
|
|||
#include "cons_data.fh"
|
||||
c
|
||||
integer n,indx(*)
|
||||
double precision coef(*),K0,gcut0,gamma0
|
||||
double precision coef(*),K0,gcut0,gamma0,c0
|
||||
double precision rion(3,*)
|
||||
double precision gamm,energy
|
||||
c
|
||||
|
|
@ -1315,15 +1339,15 @@ c
|
|||
dist = dsqrt(x*x + y*y + z*z)
|
||||
gamm = gamm + coef(i)*dist
|
||||
end do
|
||||
energy = K0*exp(-((gamm-gamma0)/gcut0)**2)
|
||||
energy = K0*exp(-((gamm-gamma0)/gcut0)**c0)
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine cons_pbondings_force(n,indx,coef,K0,gcut0,gamma0,rion,
|
||||
> gamm,energy,f)
|
||||
subroutine cons_pbondings_force(n,indx,coef,K0,gcut0,gamma0,c0,
|
||||
> rion,gamm,energy,f)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
|
|
@ -1332,7 +1356,7 @@ c
|
|||
#include "cons_data.fh"
|
||||
c
|
||||
integer n,indx(*)
|
||||
double precision coef(*),K0,gcut0,gamma0
|
||||
double precision coef(*),K0,gcut0,gamma0,c0
|
||||
double precision rion(3,*)
|
||||
double precision gamm,energy,f(3,*)
|
||||
c
|
||||
|
|
@ -1350,10 +1374,8 @@ c
|
|||
dist = dsqrt(x*x + y*y + z*z)
|
||||
gamm = gamm + coef(i)*dist
|
||||
end do
|
||||
energy = K0*exp(-((gamm-gamma0)/gcut0)**2)
|
||||
!denergy = 2.0d0*K0*(gamm-gamma0)
|
||||
denergy = -2.0d0*K0*(gamm-gamma0)*exp(-((gamm-gamma0)/gcut0)**2)
|
||||
> / gcut0**2
|
||||
energy = K0*exp(-((gamm-gamma0)/gcut0)**c0)
|
||||
denergy = -(c0*energy*(gamm-gamma0)**(c0-1.0d0)) / (gcut0**c0)
|
||||
|
||||
do i=1,n
|
||||
i0 = indx(2*i-1)
|
||||
|
|
|
|||
|
|
@ -370,7 +370,7 @@ c
|
|||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_indx"
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
subroutine cons_pbondings_coef_tag(cons_name,tag)
|
||||
|
|
@ -385,7 +385,7 @@ c
|
|||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_coef"
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
subroutine cons_pbondings_k0_tag(cons_name,tag)
|
||||
|
|
@ -400,7 +400,7 @@ c
|
|||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_k0"
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
subroutine cons_pbondings_gcut0_tag(cons_name,tag)
|
||||
|
|
@ -415,7 +415,7 @@ c
|
|||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_gcut0"
|
||||
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
|
@ -431,7 +431,22 @@ c
|
|||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_gamma0"
|
||||
return
|
||||
end
|
||||
|
||||
subroutine cons_pbondings_c0_tag(cons_name,tag)
|
||||
implicit none
|
||||
#include "util.fh"
|
||||
#include "inp.fh"
|
||||
|
||||
character*(*) tag
|
||||
c
|
||||
integer n
|
||||
character*255 cons_name
|
||||
c
|
||||
n=inp_strlen(cons_name)
|
||||
tag="cons:"//cons_name(1:n)//":pbondings_c0"
|
||||
return
|
||||
end
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue