adding Debey screening and old torsion parm

[unres4.git] / source / unres / energy.F90

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index 0e98d45..1d24115 100644 (file)
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -25514,7 +25514,7 @@
       use calc_data
       use comm_momo
        real (kind=8) ::  facd3, facd4, federmaus, adler,&
-         Ecl,Egb,Epol,Fisocav,Elj,Fgb
+         Ecl,Egb,Epol,Fisocav,Elj,Fgb,debkap
 !       integer :: k
 !c! Epol and Gpol analytical parameters
        alphapol1 = alphapol(itypi,itypj)
@@ -25559,10 +25559,15 @@
        dGCLdOM12 = 0.0d0
        ee0 = dexp(-( Rhead_sq ) / (4.0d0 * a12sq))
        Fgb = sqrt( ( Rhead_sq ) + a12sq * ee0)
-       Egb = -(332.0d0 * Qij * eps_inout_fac) / Fgb
+       debkap=debaykap(itypi,itypj)
+       Egb = -(332.0d0 * Qij *&
+        (1.0/eps_in-dexp(-debkap*Fgb)/eps_out)) / Fgb
 !       print *,"EGB WTF",Qij,eps_inout_fac,Fgb,itypi,itypj,eps_in,eps_out
 !c! Derivative of Egb is Ggb...
-       dGGBdFGB = -(-332.0d0 * Qij * eps_inout_fac) / (Fgb * Fgb)
+       dGGBdFGB = -(-332.0d0 * Qij * &
+       (1.0/eps_in-dexp(-debkap*Fgb)/eps_out))/(Fgb*Fgb)&
+       -(332.0d0 * Qij *&
+        (dexp(-debkap*Fgb)*debkap/eps_out))/ Fgb
        dFGBdR = ( Rhead * ( 2.0d0 - (0.5d0 * ee0) ) )/ ( 2.0d0 * Fgb )
        dGGBdR = dGGBdFGB * dFGBdR
 !c!-------------------------------------------------------------------