Commit Adam 6/29/2014

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

diff --git a/source/wham/src-M-NEWCORR/energy_p_new.F b/source/wham/src-M-NEWCORR/energy_p_new.F
index 0dddd7f..917b6a6 100644 (file)
--- a/source/wham/src-M-NEWCORR/energy_p_new.F
+++ b/source/wham/src-M-NEWCORR/energy_p_new.F
@@ -1718,8 +1718,8 @@ c       write (iout,*) 'theta=', theta(i-1)
           iti1=ntortyp+1
         endif
 cd        write (iout,*) '*******i',i,' iti1',iti
-cd        write (iout,*) 'b1',b1(:,iti)
-cd        write (iout,*) 'b2',b2(:,iti)
+cd        write (iout,*) 'b1',b1(:,i-2)
+cd        write (iout,*) 'b2',b2(:,i-2)
 cd        write (iout,*) 'Ug',Ug(:,:,i-2)
 c        print *,"itilde1 i iti iti1",i,iti,iti1
         if (i .gt. iatel_s+2) then
@@ -4362,7 +4362,6 @@ C-----------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
-      logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
@@ -4447,10 +4446,11 @@ c------------------------------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
-      logical lprn
+      logical lprn,energy_dec
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
+      energy_dec=.false.
       etors=0.0D0
       do i=iphi_start,iphi_end
         if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
@@ -4465,6 +4465,7 @@ c      lprn=.true.
         itori1=itortyp(itype(i-1))
         phii=phi(i)
         gloci=0.0D0
+        etors_ii=0.0d0
 C Regular cosine and sine terms
         do j=1,nterm(itori,itori1,iblock)
           v1ij=v1(j,itori,itori1,iblock)
@@ -4472,6 +4473,8 @@ C Regular cosine and sine terms
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
           etors=etors+v1ij*cosphi+v2ij*sinphi
+          if (energy_dec) etors_ii=etors_ii+
+     &      v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
 C Lorentz terms
@@ -4488,16 +4491,19 @@ C
           pom=vl2ij*cosphi+vl3ij*sinphi
           pom1=1.0d0/(pom*pom+1.0d0)
           etors=etors+vl1ij*pom1
-c          if (energy_dec) etors_ii=etors_ii+
-c     &                vl1ij*pom1
+          if (energy_dec) etors_ii=etors_ii+
+     &                vl1ij*pom1
           pom=-pom*pom1*pom1
           gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
         enddo
 C Subtract the constant term
         etors=etors-v0(itori,itori1,iblock)
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &       'etor',i,etors_ii
         if (lprn)
-     &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
+     &  write (iout,'(2(a3,2x,i3,2x),2i3,f10.2,6f8.3/36x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
+     &  rad2deg*phii,
      &  (v1(j,itori,itori1,1),j=1,6),(v2(j,itori,itori1,1),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
@@ -4550,8 +4556,8 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.TORCNSTR'
       logical lprn
 C Set lprn=.true. for debugging
-      lprn=.false.
-c     lprn=.true.
+c      lprn=.false.
+      lprn=.true.
       etors_d=0.0D0
       do i=iphi_start,iphi_end-1
         if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1