Added wham, cluster for NEWCORR

diff --git a/source/unres/src_MD-M-newcorr/Makefile_MPICH_ifort b/source/unres/src_MD-M-newcorr/Makefile_MPICH_ifort
index 708cf8f..fbc0a98 100644 (file)
--- a/source/unres/src_MD-M-newcorr/Makefile_MPICH_ifort
+++ b/source/unres/src_MD-M-newcorr/Makefile_MPICH_ifort
@@ -70,7 +70,7 @@ E0LL2Y: ${object} xdrf/libxdrf.a
        ${FC} ${OPT} ${object} cinfo.o ${LIBS}  -o ${BIN}
 
 E0LL2YT: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI \
-       -DSPLITELE -DLANG0 -DNEWCORR
+       -DSPLITELE -DLANG0 -DNEWCORR #-DMUOUT
 E0LL2YT: BIN = ../../../bin/unres/MD/unres_mchain_ifort_MPICH_E0LL2Y-NEWC.exe
 E0LL2YT: ${object} xdrf/libxdrf.a
        cc -o compinfo compinfo.c

diff --git a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F
index 42005e1..97c4144 100644 (file)
--- a/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M-newcorr/energy_p_new_barrier.F
@@ -2271,27 +2271,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 +2308,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 +2323,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 +2470,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)

diff --git a/source/unres/src_MD-M-newcorr/parmread.F b/source/unres/src_MD-M-newcorr/parmread.F
index 03efdb2..54345a5 100644 (file)
--- a/source/unres/src_MD-M-newcorr/parmread.F
+++ b/source/unres/src_MD-M-newcorr/parmread.F
@@ -352,6 +352,7 @@ C Control printout of the coefficients of virtual-bond-angle potentials
 C
       if (lprint) then
         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
+        do iblock=1,2
         do i=1,nthetyp+1
           do j=1,nthetyp+1
             do k=1,nthetyp+1
@@ -385,6 +386,7 @@ C
             enddo
           enddo
         enddo
+        enddo ! iblock
       enddo
       call flush(iout)
       endif
@@ -652,6 +654,7 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
       close (itorp)
       if (lprint) then
         write (iout,'(/a/)') 'Torsional constants:'
+        do iblock=1,2
         do i=1,ntortyp
           do j=1,ntortyp
             write (iout,*) 'ityp',i,' jtyp',j
@@ -667,6 +670,7 @@ c      &v2(k,-i,-j,iblock),v2(k,i,j,iblock)
             enddo
           enddo
         enddo
+        enddo
       endif
 
 C
@@ -887,7 +891,8 @@ cc maxinter is maximum interaction sites
             write (iout,*) 'ityp',i,' jtyp',j
             write (iout,*) 'Fourier constants'
             do k=1,nterm_sccor(i,j)
-      write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
+      write (iout,'(2(1pe15.5))') (v1sccor(k,l,i,j),v2sccor(k,l,i,j),
+     &   l=1,maxinter)
             enddo
             write (iout,*) 'Lorenz constants'
             do k=1,nlor_sccor(i,j)
@@ -1010,10 +1015,14 @@ c      lprint=.true.
       if (lprint) then
       do i=1,nloctyp
         write (iout,*) 'Type',i
-        write (iout,*) 'B1'
-        write(iout,*) B1(1,i),B1(2,i)
-        write (iout,*) 'B2'
-        write(iout,*) B2(1,i),B2(2,i)
+        write (iout,*) 'BNEW1(1)'
+        write (iout,*) (BNEW1(ii,1,i),ii=1,3)
+        write (iout,*) 'BNEW1(2)'
+        write (iout,*) BNEW1(1,2,i)
+        write (iout,*) 'BNEW2(1)'
+        write (iout,*) (BNEW2(ii,1,i),ii=1,3)
+        write (iout,*) 'BNEW2(2)'
+        write (iout,*) BNEW2(1,2,i)
         write (iout,*) 'CC'
         do j=1,2
           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
@@ -1075,7 +1084,7 @@ C
      & ', exponents are ',expon,2*expon 
       goto (10,20,30,30,40) ipot
 C----------------------- LJ potential ---------------------------------
-   10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &   (sigma0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJ potential:'
@@ -1087,7 +1096,7 @@ C----------------------- LJ potential ---------------------------------
       endif
       goto 50
 C----------------------- LJK potential --------------------------------
-   20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
       if (lprint) then
        write (iout,'(/a/)') 'Parameters of the LJK potential:'
@@ -1101,12 +1110,12 @@ C----------------------- LJK potential --------------------------------
       goto 50
 C---------------------- GB or BP potential -----------------------------
    30 do i=1,ntyp
-       read (isidep,*,end=116,err=116)(eps(i,j),j=i,ntyp)
+       read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
       enddo
-      read (isidep,*,end=116,err=116)(sigma0(i),i=1,ntyp)
-      read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
-      read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
-      read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
+      read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+      read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
+      read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
+      read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
 C For the GB potential convert sigma'**2 into chi'
       if (ipot.eq.4) then
        do i=1,ntyp
@@ -1125,7 +1134,7 @@ C For the GB potential convert sigma'**2 into chi'
       endif
       goto 50
 C--------------------- GBV potential -----------------------------------
-   40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
+   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
       if (lprint) then