From 65810a6726927aa0c848d2b68cecdae5817a40d0 Mon Sep 17 00:00:00 2001 From: Eric Bylaska Date: Sun, 7 Jun 2020 16:00:30 -0700 Subject: [PATCH] added c0 parameter to penalty function. Was Vpenalty = K0*exp(-((gamm-gamm0)/gcut0)**2) now is K0*exp(-((gamm-gamm0)/gcut0)**c0)...EJB --- src/cons/cons.F | 13 +++++--- src/cons/cons.fh | 5 +++ src/cons/cons_data.F | 42 ++++++++++++++++++++++-- src/cons/cons_data.fh | 3 +- src/cons/cons_input.F | 17 ++++++---- src/cons/cons_rtdb.F | 47 +++++++++++++++++++++++---- src/cons/cons_springs.F | 72 +++++++++++++++++++++++++++-------------- src/cons/cons_utils.F | 23 ++++++++++--- 8 files changed, 171 insertions(+), 51 deletions(-) diff --git a/src/cons/cons.F b/src/cons/cons.F index 121bb7abe5..a5a1818904 100644 --- a/src/cons/cons.F +++ b/src/cons/cons.F @@ -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)') diff --git a/src/cons/cons.fh b/src/cons/cons.fh index 196eeb8ccb..6b98fba207 100644 --- a/src/cons/cons.fh +++ b/src/cons/cons.fh @@ -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 diff --git a/src/cons/cons_data.F b/src/cons/cons_data.F index a312396062..2dd80cfdb6 100644 --- a/src/cons/cons_data.F +++ b/src/cons/cons_data.F @@ -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 + diff --git a/src/cons/cons_data.fh b/src/cons/cons_data.fh index b7eb5f7778..073fe92dba 100644 --- a/src/cons/cons_data.fh +++ b/src/cons/cons_data.fh @@ -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, diff --git a/src/cons/cons_input.F b/src/cons/cons_input.F index 7e12de4d92..b9879e0e3c 100644 --- a/src/cons/cons_input.F +++ b/src/cons/cons_input.F @@ -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 diff --git a/src/cons/cons_rtdb.F b/src/cons/cons_rtdb.F index 5c68251d7b..7115e59468 100644 --- a/src/cons/cons_rtdb.F +++ b/src/cons/cons_rtdb.F @@ -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) diff --git a/src/cons/cons_springs.F b/src/cons/cons_springs.F index d0b96a1d4f..14f40930d4 100644 --- a/src/cons/cons_springs.F +++ b/src/cons/cons_springs.F @@ -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) diff --git a/src/cons/cons_utils.F b/src/cons/cons_utils.F index 231ebf4ff5..19a5b8b2f3 100644 --- a/src/cons/cons_utils.F +++ b/src/cons/cons_utils.F @@ -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