debugging 5Dia
[unres4.git] / source / unres / energy.f90
index 4678cb5..0dbcd9c 100644 (file)
       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
         difi=phii-phi0(i)
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         endif
 !        write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
 !     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
         difi=pinorm(phii-phi0(i))
         if (difi.gt.drange(i)) then
           difi=difi-drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else if (difi.lt.-drange(i)) then
           difi=difi+drange(i)
-          edihcnstr=edihcnstr+0.25d0*ftors*difi**4
-          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+          edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4
+          gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3
         else
           difi=0.0
         endif
       write(iout,*) 'Calling CHECK_ECARTINT.'
       nf=0
       icall=0
-      write (iout,*) "Before geom_to_var"
       call geom_to_var(nvar,x)
-      write (iout,*) "after geom_to_var"
       write (iout,*) "split_ene ",split_ene
       call flush(iout)
       if (.not.split_ene) then
-        write(iout,*) 'Calling CHECK_ECARTINT if'
         call etotal(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
         etot=energia(0)
-        write (iout,*) "etot",etot
-        call flush(iout)
-!el        call enerprint(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
         do i=1,nres
           write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
         do j=1,3
           grad_s(j,0)=gcart(j,0)
         enddo
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
         do i=1,nres
           do j=1,3
             grad_s(j,i)=gcart(j,i)
           enddo
         enddo
       else
-write(iout,*) 'Calling CHECK_ECARTIN else.'
 !- split gradient check
         call zerograd
         call etotal_long(energia)
 !el        call enerprint(energia)
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
-        write (iout,*) "longrange grad"
         do i=1,nres
           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
           (gxcart(j,i),j=1,3)
@@ -11858,14 +11837,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         call zerograd
         call etotal_short(energia)
         call enerprint(energia)
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
-        write (iout,*) "shortrange grad"
         do i=1,nres
           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
           (gxcart(j,i),j=1,3)
@@ -12048,12 +12021,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         call etotal(energia)
         etot=energia(0)
 !el        call enerprint(energia)
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
         do i=1,nres
           write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
@@ -12073,14 +12041,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         call zerograd
         call etotal_long(energia)
 !el        call enerprint(energia)
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
-        write (iout,*) "longrange grad"
         do i=1,nres
           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
           (gxcart(j,i),j=1,3)
@@ -12097,14 +12059,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         call zerograd
         call etotal_short(energia)
 !el        call enerprint(energia)
-        call flush(iout)
-        write (iout,*) "enter cartgrad"
-        call flush(iout)
         call cartgrad
-        write (iout,*) "exit cartgrad"
-        call flush(iout)
         icall =1
-        write (iout,*) "shortrange grad"
         do i=1,nres
           write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
           (gxcart(j,i),j=1,3)
@@ -16230,7 +16186,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 ! This subrouting calculates total Cartesian coordinate gradient. 
 ! The subroutine chainbuild_cart and energy MUST be called beforehand.
 !
-!el#define DEBUG
+#define DEBUG
 #ifdef TIMING
       time00=MPI_Wtime()
 #endif
@@ -16332,7 +16288,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 #ifdef TIMING
             time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 #endif
-      !el#undef DEBUG
+#undef DEBUG
             return
             end subroutine cartgrad
       !-----------------------------------------------------------------------------
@@ -20419,7 +20375,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 !c
 !c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 !c
-      print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
+!      print *,"iatel_s_nucl,iatel_e_nucl",iatel_s_nucl,iatel_e_nucl
       do i=iatel_s_nucl,iatel_e_nucl
         if (itype(i,2).eq.ntyp1_molec(2) .or. itype(i+1,2).eq.ntyp1_molec(2)) cycle
         dxi=dc(1,i)
@@ -20806,8 +20762,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, &
@@ -21764,14 +21720,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)
@@ -21790,6 +21749,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
@@ -21847,6 +21807,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