MPICH2 and energy fixes in WHAM-M and Cluster-M.

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

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index c588bee..4814c1d 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -124,7 +124,6 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
 #endif
       energia(0)=etot
       energia(1)=evdw
-c      call enerprint(energia(0),frac)
 #ifdef SCP14
       energia(2)=evdw2-evdw2_14
       energia(17)=evdw2_14
@@ -229,7 +228,8 @@ C
      &   +wturn3*fact(2)*gel_loc_turn3(i)
      &   +wturn6*fact(5)*gel_loc_turn6(i)
      &   +wel_loc*fact(2)*gel_loc_loc(i)
-     &   +wsccor*fact(1)*gsccor_loc(i)
+c     &   +wsccor*fact(1)*gsccor_loc(i)
+c ROZNICA Z WHAMem
       enddo
       endif
       return
@@ -360,6 +360,14 @@ C
       integer icant
       external icant
 cd    print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+c ROZNICA DODANE Z WHAM
+c      do i=1,210
+c        do j=1,2
+c          eneps_temp(j,i)=0.0d0
+c        enddo
+c      enddo
+cROZNICA
+
       evdw=0.0D0
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
@@ -393,6 +401,11 @@ c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             e2=fac*bb(itypi,itypj)
             evdwij=e1+e2
             ij=icant(itypi,itypj)
+c ROZNICA z WHAM
+c            eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
+c            eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij
+c
+
 cd          sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
 cd          epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
 cd          write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
@@ -3122,8 +3135,8 @@ C
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
       double precision y(2),z(2)
       delta=0.02d0*pi
-      time11=dexp(-2*time)
-      time12=1.0d0
+c      time11=dexp(-2*time)
+c      time12=1.0d0
       etheta=0.0D0
 c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
@@ -3148,8 +3161,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
-          icrc=0
-          call proc_proc(phii,icrc)
+c          icrc=0
+c          call proc_proc(phii,icrc)
           if (icrc.eq.1) phii=150.0
 #else
           phii=phi(i)
@@ -3163,8 +3176,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.lt.nres .and. itype(i).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
-          icrc=0
-          call proc_proc(phii1,icrc)
+c          icrc=0
+c          call proc_proc(phii1,icrc)
           if (icrc.eq.1) phii1=150.0
           phii1=pinorm(phii1)
           z(1)=cos(phii1)
@@ -3232,7 +3245,7 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
- 1215   continue
+c 1215   continue
       enddo
 C Ufff.... We've done all this!!! 
       return
@@ -3375,7 +3388,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+c        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.(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
@@ -3388,7 +3403,7 @@ CC Ta zmina jest niewlasciwa
           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(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3402,13 +3417,14 @@ CC Ta zmina jest niewlasciwa
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+c          ityp1=nthetyp+1
           do k=1,nsingle
+            ityp1=ithetyp((itype(i-2)))
             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
@@ -3423,7 +3439,8 @@ CC Ta zmina jest niewlasciwa
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+c          ityp3=nthetyp+1
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -3540,7 +3557,8 @@ c        call flush(iout)
         etheta=etheta+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
+c        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
       return
       end
@@ -3931,7 +3949,8 @@ c        write (2,*) "xx",xx," yy",yy," zz",zz
 Cc diagnostics - remove later
         xx1 = dcos(alph(2))
         yy1 = dsin(alph(2))*dcos(omeg(2))
-        zz1 = -dsin(alph(2))*dsin(omeg(2))
+c        zz1 = -dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0d0,itype(i))*dsin(alph(2))*dsin(omeg(2))
         write(2,'(3f8.1,3f9.3,1x,3f9.3)') 
      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
      &    xx1,yy1,zz1
@@ -4524,7 +4543,7 @@ c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
       do i=itau_start,itau_end
-        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
         esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
@@ -5233,7 +5252,11 @@ C---------------------------------------------------------------------------
      &  auxmat(2,2)
       iti1 = itortyp(itype(i+1))
       if (j.lt.nres-1) then
-        itj1 = itortyp(itype(j+1))
+        if (itype(j).le.ntyp) then
+          itj1 = itortyp(itype(j+1))
+        else
+          itj1=ntortyp+1
+        endif
       else
         itj1=ntortyp+1
       endif
@@ -5321,14 +5344,16 @@ cd      if (i.ne.2 .or. j.ne.4 .or. k.ne.3 .or. l.ne.5) return
       enddo 
       if (l.eq.j+1) then
 C parallel orientation of the two CA-CA-CA frames.
-        if (i.gt.1) then
+c        if (i.gt.1) then
+        if (i.gt.1 .and. itype(i).le.ntyp) then
           iti=itortyp(itype(i))
         else
           iti=ntortyp+1
         endif
         itk1=itortyp(itype(k+1))
         itj=itortyp(itype(j))
-        if (l.lt.nres-1) then
+c        if (l.lt.nres-1) then
+        if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then
           itl1=itortyp(itype(l+1))
         else
           itl1=ntortyp+1
@@ -5474,7 +5499,8 @@ C Calculate the Cartesian derivatives of the vectors.
 C End vectors
       else
 C Antiparallel orientation of the two CA-CA-CA frames.
-        if (i.gt.1) then
+c        if (i.gt.1) then
+        if (i.gt.1 .and. itype(i).le.ntyp) then
           iti=itortyp(itype(i))
         else
           iti=ntortyp+1
@@ -5482,7 +5508,8 @@ C Antiparallel orientation of the two CA-CA-CA frames.
         itk1=itortyp(itype(k+1))
         itl=itortyp(itype(l))
         itj=itortyp(itype(j))
-        if (j.lt.nres-1) then
+c        if (j.lt.nres-1) then
+        if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
           itj1=itortyp(itype(j+1))
         else 
           itj1=ntortyp+1
@@ -6640,14 +6667,16 @@ C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 C           energy moment and not to the cluster cumulant.
       iti=itortyp(itype(i))
-      if (j.lt.nres-1) then
+c      if (j.lt.nres-1) then
+      if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
         itj1=itortyp(itype(j+1))
       else
         itj1=ntortyp+1
       endif
       itk=itortyp(itype(k))
       itk1=itortyp(itype(k+1))
-      if (l.lt.nres-1) then
+c      if (l.lt.nres-1) then
+      if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then
         itl1=itortyp(itype(l+1))
       else
         itl1=ntortyp+1
@@ -6760,13 +6789,15 @@ C           energy moment and not to the cluster cumulant.
 cd      write (2,*) 'eello_graph4: wturn6',wturn6
       iti=itortyp(itype(i))
       itj=itortyp(itype(j))
-      if (j.lt.nres-1) then
+c      if (j.lt.nres-1) then
+      if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then
         itj1=itortyp(itype(j+1))
       else
         itj1=ntortyp+1
       endif
       itk=itortyp(itype(k))
-      if (k.lt.nres-1) then
+c      if (k.lt.nres-1) then
+      if (k.lt.nres-1 .and. itype(k+1).le.ntyp) then
         itk1=itortyp(itype(k+1))
       else
         itk1=ntortyp+1