energy match prerelease from old unres
[unres4.git] / source / unres / energy.f90
index 6bb718a..4bb6804 100644 (file)
 !d     &      xmedi,ymedi,zmedi,xj,yj,zj
 
           if (energy_dec) then 
-              write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &
-                  'evdw1',i,j,evdwij,&
-                  iteli,itelj,aaa,evdw1
+!              write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &
+!                  'evdw1',i,j,evdwij,&
+!                  iteli,itelj,aaa,evdw1
+              write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
               write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
           endif
 !
           endif
           evdwij=e1+e2
           evdw2=evdw2+evdwij
-          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
-             'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
-            bad(itypj,iteli)
+!          if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
+!             'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
+          if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+             'evdw2',i,j,evdwij
 !
 ! Calculate contributions to the gradient in the virtual-bond and SC vectors.
 !
              "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
         else
         diff = vbld(i)-vbldp0
-        if (energy_dec) write (iout,*) &
+        if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
            "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
       etheta=0.0D0
       do i=ithet_start,ithet_end
         if (itype(i-1).eq.ntyp1) cycle
+        if (itype(i-2).eq.ntyp1.or.itype(i).eq.ntyp1) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
          phii1*rad2deg,ethetai
 !        lprn1=.false.
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+                                    'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
         gloc(nphi+i-2,icg)=wang*dethetai
             gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
           enddo
         endif
-        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+        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/)') &
 !      include 'COMMON.IOUNITS'
 !      include 'COMMON.FFIELD'
 !      include 'COMMON.TORCNSTR'
-      real(kind=8) :: etors_d
+      real(kind=8) :: etors_d,etors_d_ii
       logical :: lprn
 !el local variables
       integer :: i,j,k,l,itori,itori1,itori2,iblock
       etors_d=0.0D0
 !      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
+        etors_d_ii=0.0D0
         if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
             .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         itori=itortyp(itype(i-2))
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ &
            v2cij*cosphi2+v2sij*sinphi2
+          if (energy_dec) etors_d_ii=etors_d_ii+ &
+           v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ &
               v1sdij*sinphi1p2+v2sdij*sinphi1m2
+            if (energy_dec) etors_d_ii=etors_d_ii+ &
+              v1cdij*cosphi1p2+v2cdij*cosphi1m2+ &
+              v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 &
               -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 &
               -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+                            'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
       enddo
 !      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         phii=phi(i)
         do intertyp=1,3 !intertyp
+         esccor_ii=0.0D0
 !c Added 09 May 2012 (Adasko)
 !c  Intertyp means interaction type of backbone mainchain correlation: 
 !   1 = SC...Ca...Ca...Ca
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
+          if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+        if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') &
+                                'esccor',i,intertyp,esccor_ii
 !      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn) &