changes in wham and unres

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

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index b867e3f..7ab625f 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -604,6 +604,7 @@
        etube=0.0d0
       endif
 !--------------------------------------------------------
+!      write (iout,*) "NRES_MOLEC(2),",nres_molec(2)
 !      print *,"before",ees,evdw1,ecorr
       if (nres_molec(2).gt.0) then
       call ebond_nucl(estr_nucl)
@@ -614,6 +615,13 @@
       call epsb(evdwpsb,eelpsb)
       call esb(esbloc)
       call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+      else
+       etors_nucl=0.0d0
+       estr_nucl=0.0d0
+       ebe_nucl=0.0d0
+       evdwsb=0.0d0
+       eelsb=0.0d0
+       esbloc=0.0d0
       endif
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
@@ -20799,8 +20807,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
             evdwij=evdwij*eps2rt
             evdwsb=evdwsb+evdwij
             if (lprn) then
-            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+            sigm=dabs(aa_nucl(itypi,itypj)/bb_nucl(itypi,itypj))**(1.0D0/6.0D0)
+            epsi=bb_nucl(itypi,itypj)**2/aa_nucl(itypi,itypj)
             write (iout,'(2(a3,i3,2x),17(0pf7.3))') &
              restyp(itypi,2),i,restyp(itypj,2),j, &
              epsi,sigm,chi1,chi2,chip1,chip2, &
@@ -21757,14 +21765,17 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         r012 = r06**2
         k0 = 332.0*(2.0*2.0)/80.0
         itmp=0
+        
         do i=1,4
         itmp=itmp+nres_molec(i)
         enddo
+        write(iout,*) "itmp",itmp
         do i=itmp+1,itmp+nres_molec(5)-1
        
         xi=c(1,i)
         yi=c(2,i)
         zi=c(3,i)
+         
           xi=mod(xi,boxxsize)
           if (xi.lt.0) xi=xi+boxxsize
           yi=mod(yi,boxysize)
@@ -21783,6 +21794,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           if (yj.lt.0) yj=yj+boxysize
           zj=dmod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
+          write(iout,*) c(1,i),xi,xj,"xy",boxxsize
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
       yj_safe=yj
@@ -21840,6 +21852,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
         enddo
 
+        write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
         ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
        enddo
        enddo
@@ -23370,7 +23383,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       real(kind=8) :: evdw,sig0ij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
                     dist_temp, dist_init,aa,bb,ssgradlipi,ssgradlipj, &
-                    sslipi,sslipj,faclip
+                    sslipi,sslipj,faclip,alpha_sco
       integer :: ii
       real(kind=8) :: fracinbuf
        real (kind=8) :: escpho
@@ -23426,6 +23439,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           do xshift=-1,1
           do yshift=-1,1
           do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
           yj=yj_safe+yshift*boxysize
           zj=zj_safe+zshift*boxzsize
           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
@@ -23652,12 +23666,15 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !c!      write (*,*) "Gvdwc(",k,",",i,")=", gvdwc(k,i)
 !c!      write (*,*) "Gvdwc(",k,",",j,")=", gvdwc(k,j)
 !      alphapol1 = alphapol_scpho(itypi)
-       if (wqq_scpho(itypi).gt.0.0) then
+       if (wqq_scpho(itypi).ne.0.0) then
        Qij=wqq_scpho(itypi)/eps_in
+       alpha_sco=1.d0/alphi_scpho(itypi)
 !       Qij=0.0
-       Ecl = (332.0d0 * Qij) / Rhead
+       Ecl = (332.0d0 * Qij*dexp(-Rhead*alpha_sco)) / Rhead
 !c! derivative of Ecl is Gcl...
-       dGCLdR = (-332.0d0 * Qij ) / Rhead_sq
+       dGCLdR = (-332.0d0 * Qij*dexp(-Rhead*alpha_sco)*  &
+                (Rhead*alpha_sco+1) ) / Rhead_sq
+       if (energy_dec) write(iout,*) "ECL",ECL,Rhead,1.0/rij
        else if (wqdip_scpho(2,itypi).gt.0.0d0) then
        w1        = wqdip_scpho(1,itypi)
        w2        = wqdip_scpho(2,itypi)
@@ -23679,6 +23696,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        Ecl = sparrow / Rhead**2.0d0 &
            - hawk    / Rhead**4.0d0
 !c!-------------------------------------------------------------------
+       if (energy_dec) write(iout,*) "ECLdipdip",ECL,Rhead,&
+           1.0/rij,sparrow
+
 !c! derivative of ecl is Gcl
 !c! dF/dr part
        dGCLdR  = - 2.0d0 * sparrow / Rhead**3.0d0 &
@@ -23780,6 +23800,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 
        END DO
 !       if (i.eq.3) print *,i,j,evdwij,epol,Fcav,ECL
+       if (energy_dec) write (iout,'(a22,2i5,4f8.3,f16.3)'), &
+        "escpho:evdw,pol,cav,CL",i,j,evdwij,epol,Fcav,ECL,escpho
        escpho=escpho+evdwij+epol+Fcav+ECL
        call sc_grad_scpho
          enddo
@@ -23927,6 +23949,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           do xshift=-1,1
           do yshift=-1,1
           do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
           yj=yj_safe+yshift*boxysize
           zj=zj_safe+zshift*boxzsize
           dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2