Adam's cluster & unres corrections

[unres.git] / source / unres / src-HCD-5D / energy_p_new_barrier.F

diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F
index ddc2a4d..2c701ca 100644 (file)
--- a/source/unres/src-HCD-5D/energy_p_new_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F
@@ -118,9 +118,17 @@ c        call chainbuild_cart
 c      write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
       if (mod(itime_mat,imatupdate).eq.0) then
         call make_SCp_inter_list
+c        write (iout,*) "Finished make_SCp_inter_list"
+c        call flush(iout)
         call make_SCSC_inter_list
+c        write (iout,*) "Finished make_SCSC_inter_list"
+c        call flush(iout)
         call make_pp_inter_list
+c        write (iout,*) "Finished make_pp_inter_list"
+c        call flush(iout)
         call make_pp_vdw_inter_list
+c        write (iout,*) "Finished make_pp_vdw_inter_list"
+c        call flush(iout)
       endif
 c      print *,'Processor',myrank,' calling etotal ipot=',ipot
 c      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
@@ -366,7 +374,17 @@ c         call flush(iout)
 c         write (iout,*) "MULTIBODY_HB ecorr",ecorr,ecorr5,ecorr6,n_corr,
 c     &     n_corr1
 c         call flush(iout)
+      else
+         ecorr=0.0d0
+         ecorr5=0.0d0
+         ecorr6=0.0d0
+         eturn6=0.0d0
       endif
+#else
+      ecorr=0.0d0
+      ecorr5=0.0d0
+      ecorr6=0.0d0
+      eturn6=0.0d0
 #endif
 c      print *,"Processor",myrank," computed Ucorr"
 c      write (iout,*) "nsaxs",nsaxs," saxs_mode",saxs_mode
@@ -409,15 +427,17 @@ C      print *,"za lipidami"
         call AFMforce(Eafmforce)
       else if (selfguide.gt.0) then
         call AFMvel(Eafmforce)
+      else 
+        Eafmforce=0.0d0
       endif
       if (TUBElog.eq.1) then
 C      print *,"just before call"
         call calctube(Etube)
-       elseif (TUBElog.eq.2) then
+      elseif (TUBElog.eq.2) then
         call calctube2(Etube)
-       else
-       Etube=0.0d0
-       endif
+      else
+        Etube=0.0d0
+      endif
 
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
@@ -2053,6 +2073,9 @@ C derivatives.
               sigsq=1.0D0/sigsq
               sig=sig0ij*dsqrt(sigsq)
               rij_shift=1.0D0/rij-sig+sig0ij
+c              if (energy_dec)
+c     &        write (iout,*) "rij",1.0d0/rij," rij_shift",rij_shift,
+c     &       " sig",sig," sig0ij",sig0ij
 c for diagnostics; uncomment
 c            rij_shift=1.2*sig0ij
 C I hate to put IF's in the loops, but here don't have another choice!!!!
@@ -2061,7 +2084,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!!
 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 cd     &        restyp(itypi),i,restyp(itypj),j,
 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
-                return
+c                return
               endif
               sigder=-sig*sigsq
 c---------------------------------------------------------------
@@ -2093,7 +2116,7 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &           evdwij
               endif
 
-              if (energy_dec) write (iout,'(a,2i5,3f10.5)') 
+              if (energy_dec) write (iout,'(a,2i5,2f10.5,e15.5)') 
      &                    'r sss evdw',i,j,1.0d0/rij,sss,evdwij
 
 C Calculate gradient components.
@@ -2285,6 +2308,12 @@ C Calculate gradient components.
            fac=rij*fac-2*expon*rrij*e_augm
            fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
 C Calculate the radial part of the gradient
+           gg_lipi(3)=eps1*(eps2rt*eps2rt)
+     &       *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+     &        (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+     &       +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+           gg_lipj(3)=ssgradlipj*gg_lipi(3)
+           gg_lipi(3)=gg_lipi(3)*ssgradlipi
            gg(1)=xj*fac
            gg(2)=yj*fac
            gg(3)=zj*fac
@@ -7691,16 +7720,20 @@ c
          enddo
          
 c         min_odl=minval(distancek)
-         do kk=1,constr_homology
-          if(l_homo(kk,ii)) then 
-            min_odl=distancek(kk)
-            exit
-          endif
-         enddo
-         do kk=1,constr_homology
-          if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) 
+         if (nexl.gt.0) then
+           min_odl=0.0d0
+         else
+           do kk=1,constr_homology
+            if(l_homo(kk,ii)) then 
+              min_odl=distancek(kk)
+              exit
+            endif
+           enddo
+           do kk=1,constr_homology
+            if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) 
      &              min_odl=distancek(kk)
-         enddo
+           enddo
+         endif
 
 c        write (iout,* )"min_odl",min_odl
 #ifdef DEBUG