Commit Adam 6/29/2014
[unres.git] / source / unres / src_MD-M-newcorr / energy_p_new_barrier.F
index 42005e1..c9140a6 100644 (file)
@@ -710,7 +710,6 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
-#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc before reduce"
       do i=1,nres
@@ -719,7 +718,6 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
         do i=1,nres
          do j=1,3
           gloc_scbuf(j,i)=gloc_sc(j,i,icg)
@@ -739,7 +737,6 @@ c      enddo
         call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
-#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc after reduce"
       do i=1,nres
@@ -748,7 +745,6 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -2271,27 +2267,27 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
           iti1=ntortyp+1
         endif
 c        write(iout,*),i
-        b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0)
-     &           +bnew1(2,1,iti)*sin(theta(i-1))
-     &           +bnew1(3,1,iti)*cos(theta(i-1)/2.0)
-        gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0
-     &             +bnew1(2,1,iti)*cos(theta(i-1))
-     &             -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0
-c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
+        b1(1,i-2)=bnew1(1,1,iti)*dsin(theta(i-1)/2.0d0)
+     &           +bnew1(2,1,iti)*dsin(theta(i-1))
+     &           +bnew1(3,1,iti)*dcos(theta(i-1)/2.0d0)
+        gtb1(1,i-2)=bnew1(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+     &             +bnew1(2,1,iti)*dcos(theta(i-1))
+     &             -bnew1(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
+c     &           +bnew1(3,1,iti)*dsin(alpha(i))*cos(beta(i))
 c     &*(cos(theta(i)/2.0)
-        b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0)
-     &           +bnew2(2,1,iti)*sin(theta(i-1))
-     &           +bnew2(3,1,iti)*cos(theta(i-1)/2.0)
-c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
+        b2(1,i-2)=bnew2(1,1,iti)*dsin(theta(i-1)/2.0d0)
+     &           +bnew2(2,1,iti)*dsin(theta(i-1))
+     &           +bnew2(3,1,iti)*dcos(theta(i-1)/2.0d0)
+c     &           +bnew2(3,1,iti)*dsin(alpha(i))*dcos(beta(i))
 c     &*(cos(theta(i)/2.0)
-        gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0
-     &             +bnew2(2,1,iti)*cos(theta(i-1))
-     &             -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+        gtb2(1,i-2)=bnew2(1,1,iti)*dcos(theta(i-1)/2.0d0)/2.0d0
+     &             +bnew2(2,1,iti)*dcos(theta(i-1))
+     &             -bnew2(3,1,iti)*dsin(theta(i-1)/2.0d0)/2.0d0
 c        if (ggb1(1,i).eq.0.0d0) then
 c        write(iout,*) 'i=',i,ggb1(1,i),
-c     &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
-c     &bnew1(2,1,iti)*cos(theta(i)),
-c     &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
+c     &bnew1(1,1,iti)*dcos(theta(i)/2.0d0)/2.0d0,
+c     &bnew1(2,1,iti)*dcos(theta(i)),
+c     &bnew1(3,1,iti)*dsin(theta(i)/2.0d0)/2.0d0
 c        endif
         b1(2,i-2)=bnew1(1,2,iti)
         gtb1(2,i-2)=0.0
@@ -2308,8 +2304,8 @@ c        endif
 c        EE(2,2,iti)=0.0d0
 c        EE(1,2,iti)=0.5d0*eenew(1,iti)
 c        EE(2,1,iti)=0.5d0*eenew(1,iti)
-c        b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
-c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
+c        b1(2,iti)=bnew1(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
+c        b2(2,iti)=bnew2(1,2,iti)*dsin(alpha(i))*dsin(beta(i))
        b1tilde(1,i-2)=b1(1,i-2)
        b1tilde(2,i-2)=-b1(2,i-2)
        b2tilde(1,i-2)=b2(1,i-2)
@@ -2323,6 +2319,17 @@ c       write (iout,*) 'theta=', theta(i-1)
       do i=3,nres+1
 #endif
 #endif
+        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          iti = itortyp(itype(i-2))
+        else
+          iti=ntortyp+1
+        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
+          iti1 = itortyp(itype(i-1))
+        else
+          iti1=ntortyp+1
+        endif
         if (i .lt. nres+1) then
           sin1=dsin(phi(i))
           cos1=dcos(phi(i))
@@ -2459,11 +2466,12 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
           mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
 #ifdef MUOUT
-        write (iout,'(2hmu,i3,3f8.1,7f10.5)') i-2,rad2deg*theta(i-1),
+        write (iout,'(2hmu,i3,3f8.1,12f10.5)') i-2,rad2deg*theta(i-1),
      &     rad2deg*theta(i),rad2deg*phi(i),mu(1,i-2),mu(2,i-2),
      &       -b2(1,i-2),b2(2,i-2),b1(1,i-2),b1(2,i-2),
      &       dsqrt(b2(1,i-1)**2+b2(2,i-1)**2)
-     &      +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2)
+     &      +dsqrt(b1(1,i-1)**2+b1(2,i-1)**2),
+     &      ((ee(l,k,i-2),l=1,2),k=1,2),eenew(1,itortyp(iti))
 #endif
 cd        write (iout,*) 'mu ',mu(:,i-2)
 cd        write (iout,*) 'mu1',mu1(:,i-2)
@@ -5796,7 +5804,7 @@ c----------------------------------------------------------------------------
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
-c     lprn=.true.
+c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
         if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 
 C Subtract the constant term
         etors=etors-v0(itori,itori1,iblock)
           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &         'etor',i,etors_ii-v0(itori,itori1,iblock)
+     &         '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,
-     &  (v1(j,itori,itori1,iblock),j=1,6),
-     &  (v2(j,itori,itori1,iblock),j=1,6)
+     &  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*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo