Adam's corrections

[unres.git] / source / wham / src-HCD / energy_p_new.F

diff --git a/source/wham/src-HCD/energy_p_new.F b/source/wham/src-HCD/energy_p_new.F
index f4dabad..efba869 100644 (file)
--- a/source/wham/src-HCD/energy_p_new.F
+++ b/source/wham/src-HCD/energy_p_new.F
@@ -20,6 +20,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.SHIELD'
       include 'COMMON.CONTROL'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.SAXS'
       double precision fact(6)
 c      write(iout, '(a,i2)')'Calling etotal ipot=',ipot
 c      call flush(iout)
@@ -158,6 +159,7 @@ c         write (iout,*) "Calling multibody_hbond"
          call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
       endif
 #endif
+c      write (iout,*) "nsaxs",nsaxs
 c      write (iout,*) "From Esaxs: Esaxs_constr",Esaxs_constr
       if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
         call e_saxs(Esaxs_constr)
@@ -179,12 +181,16 @@ c      write(iout,*) "TEST_ENE1 constr_homology=",constr_homology
 c      write(iout,*) "TEST_ENE1 ehomology_constr=",ehomology_constr
 #ifdef DFA
 C     BARTEK for dfa test!
+      edfadis=0.0d0
       if (wdfa_dist.gt.0) call edfad(edfadis)
 c      write(iout,*)'edfad is finished!', wdfa_dist,edfadis
+      edfator=0.0d0
       if (wdfa_tor.gt.0) call edfat(edfator)
 c      write(iout,*)'edfat is finished!', wdfa_tor,edfator
+      edfanei=0.0d0
       if (wdfa_nei.gt.0) call edfan(edfanei)
 c      write(iout,*)'edfan is finished!', wdfa_nei,edfanei
+      edfabet=0.0d0
       if (wdfa_beta.gt.0) call edfab(edfabet)
 c      write(iout,*)'edfab is finished!', wdfa_beta,edfabet
 #endif
@@ -1342,8 +1348,8 @@ c#define DEBUG
 #endif
 c#undef DEBUG
 c            endif
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
-     &                        'evdw',i,j,evdwij
+            if (energy_dec) write (iout,'(a,2i5,3f10.5)')
+     &                    'r sss evdw',i,j,1.0d0/rij,sss,evdwij
             if (calc_grad) then
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
@@ -1794,13 +1800,19 @@ C to calculate the el-loc multibody terms of various order.
 C
 c      write(iout,*) 'SET_MATRICES nphi=',nphi,nres
       do i=3,nres+1
-        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        ii=ireschain(i-2)
+        if (ii.eq.0) cycle
+        innt=chain_border(1,ii)
+        inct=chain_border(2,ii)
+c        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        if (i.gt. innt+2 .and. i.lt.inct+2) then
           iti = itype2loc(itype(i-2))
         else
           iti=nloctyp
         endif
 c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
-        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+c        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+        if (i.gt. innt+1 .and. i.lt.inct+1) then
           iti1 = itype2loc(itype(i-1))
         else
           iti1=nloctyp
@@ -4598,16 +4610,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
          write (iout,*) "ij dij",i,j,dij
@@ -5165,6 +5181,10 @@ c
       estr1=0.0d0
 c      write (iout,*) "distchainmax",distchainmax
       do i=nnt+1,nct
+#ifdef FIVEDIAG
+        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+        diff = vbld(i)-vbldp0
+#else
         if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
 C          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
 C          do j=1,3
@@ -5182,13 +5202,14 @@ C         write(iout,*) i,diff
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
          endif
+#endif
           estr=estr+diff*diff
           do j=1,3
             gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
           enddo
 C        endif
-C        write (iout,'(a7,i5,4f7.3)')
-C     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
+          if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
       enddo
       estr=0.5d0*AKP*estr+estr1
 c
@@ -5200,8 +5221,9 @@ c
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
             diff=vbld(i+nres)-vbldsc0(1,iti)
-C            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-C     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec) write (iout,*) "estr sc",iti,vbld(i+nres),
+     &      vbldsc0(1,iti),diff,
+     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)