From 79213f3121536c142e03e6e2272a70381f7d8408 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Mon, 2 Nov 2015 10:26:17 +0100 Subject: [PATCH] rozgrzebany shield --- source/unres/src_MD-M/energy_p_new_barrier.F | 31 ++++++++++++++++---------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 82d6a3f..5cb9a2a 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -3708,8 +3708,8 @@ C 12/26/95 - for the evaluation of multi-body H-bonding interactions if (shield_mode.gt.0) then C fac_shield(i)=0.4 C fac_shield(j)=0.6 - el1=el1*fac_shield(i)*fac_shield(j) - el2=el2*fac_shield(i)*fac_shield(j) + el1=el1*fac_shield(i)**2*fac_shield(j)**2 + el2=el2*fac_shield(i)**2*fac_shield(j)**2 eesij=(el1+el2) ees=ees+eesij else @@ -3756,9 +3756,10 @@ C print *,i,j iresshield=shield_list(ilist,i) do k=1,3 rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i) + & *2.0 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) + & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield) C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i) @@ -3781,9 +3782,10 @@ C endif iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) + & *2.0 gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) @@ -3806,13 +3808,13 @@ C endif do k=1,3 gshieldc(k,i)=gshieldc(k,i)+ - & grad_shield(k,i)*eesij/fac_shield(i) + & grad_shield(k,i)*eesij/fac_shield(i)*2.0 gshieldc(k,j)=gshieldc(k,j)+ - & grad_shield(k,j)*eesij/fac_shield(j) + & grad_shield(k,j)*eesij/fac_shield(j)*2.0 gshieldc(k,i-1)=gshieldc(k,i-1)+ - & grad_shield(k,i)*eesij/fac_shield(i) + & grad_shield(k,i)*eesij/fac_shield(i)*2.0 gshieldc(k,j-1)=gshieldc(k,j-1)+ - & grad_shield(k,j)*eesij/fac_shield(j) + & grad_shield(k,j)*eesij/fac_shield(j)*2.0 enddo endif @@ -3931,7 +3933,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)*fac_shield(j) + & fac_shield(i)**2*fac_shield(j)**2 enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -3952,11 +3954,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)**2*fac_shield(j)**2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)**2*fac_shield(j)**2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) enddo @@ -7302,6 +7304,7 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -7720,6 +7723,7 @@ C This subroutine calculates multi-body contributions to hydrogen-bonding include 'COMMON.CONTACTS' include 'COMMON.CHAIN' include 'COMMON.CONTROL' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) integer num_cont_hb_old(maxres) logical lprn,ldone @@ -8022,6 +8026,7 @@ cd write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1) call calc_eello(i,jp,i+1,jp1,jj,kk) if (wcorr4.gt.0.0d0) & ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk) +CC & *fac_shield(i)**2*fac_shield(j)**2 if (energy_dec.and.wcorr4.gt.0.0d0) 1 write (iout,'(a6,4i5,0pf7.3)') 2 'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk) @@ -8141,6 +8146,7 @@ c------------------------------------------------------------------------------ include 'COMMON.DERIV' include 'COMMON.INTERACT' include 'COMMON.CONTACTS' + include 'COMMON.SHIELD' double precision gx(3),gx1(3) logical lprn lprn=.false. @@ -8151,7 +8157,8 @@ c------------------------------------------------------------------------------ ees0mij=ees0m(jj,i) ees0mkl=ees0m(kk,k) ekont=eij*ekl - ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl) + ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)* + & fac_shield(i)**2*fac_shield(j)**2 cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl) C Following 4 lines for diagnostics. cd ees0pkl=0.0D0 -- 1.7.9.5