MPICH2 and energy fixes in WHAM-M and Cluster-M.
authorPawel Krupa <vetinari@piasek4.chem.univ.gda.pl>
Tue, 21 Apr 2015 23:18:13 +0000 (01:18 +0200)
committerPawel Krupa <vetinari@piasek4.chem.univ.gda.pl>
Tue, 21 Apr 2015 23:18:13 +0000 (01:18 +0200)
source/cluster/wham/src-M/energy_p_new.F
source/cluster/wham/src-M/main_clust.F
source/cluster/wham/src-M/probabl.F
source/unres/src_MD-M/ssMD.F
source/wham/src-M/energy_p_new.F

index c588bee..4814c1d 100644 (file)
@@ -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
 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
index 15e0bd0..892a6c7 100644 (file)
@@ -29,12 +29,12 @@ C
       INTEGER IA(maxconf),IB(maxconf)
       INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
       INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
-      integer nn,ndis
-      real*4 DISNN
+      integer nn,ndis,scount_buf
+      real*4 DISNN, diss_buf(maxdist)
       DIMENSION NN(maxconf),DISNN(maxconf)
       LOGICAL FLAG(maxconf)
       integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
-     & it,ncon_work,ind1,kkk
+     & it,ncon_work,ind1,kkk, ijk
       double precision t1,t2,tcpu,difconf
       
       double precision varia(maxvar)
@@ -165,8 +165,16 @@ c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
       t1=tcpu()
       PRINT '(a)','End of distance computation'
 
+      scount_buf=scount(me)
+
+      do ijk=1, ndis
+      diss_buf(ijk)=diss(ijk)
+      enddo
+
+
 #ifdef MPI
-      call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
+      WRITE (iout,*) "Wchodze do call MPI_Gatherv"
+      call MPI_Gatherv(diss_buf(1),scount_buf,MPI_REAL,diss(1),
      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
       if (me.eq.master) then
 #endif
index 293fb8f..543422a 100644 (file)
       character*80 bxname
       character*2 licz1
       character*5 ctemper
-      integer ilen
+      integer ilen,ijk
       external ilen
-      real*4 Fdimless(maxconf)
-      double precision energia(0:max_ene)
+      real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
+      double precision energia(0:max_ene), totfree_buf(0:maxconf),
+     &  entfac_buf(maxconf)
       do i=1,ncon
         list_conf(i)=i
       enddo
@@ -134,6 +135,7 @@ c        call flush(iout)
           call int_from_cart1(.false.)
           call etotal(energia(0),fT)
           totfree(i)=energia(0)         
+          totfree_buf(i)=totfree(i)
 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 c          call enerprint(energia(0),fT)
@@ -208,27 +210,49 @@ c#endif
         write (iout,*) "ebe" ebe,wang
 #endif        
         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
+        Fdimless_buf(i)=Fdimless(i)
         totfree(i)=etot
+        totfree_buf(i)=totfree(i)
 #ifdef DEBUG
         write (iout,*) "fdim calc", i,ii,ib,
      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
      &   entfac(ii),Fdimless(i)
 #endif
       enddo   ! i
+
+      do ijk=1,maxconf
+      entfac_buf(ijk)=entfac(ijk)
+      Fdimless_buf(ijk)=Fdimless(ijk)
+      enddo
+      do ijk=0,maxconf
+      totfree_buf(ijk)=totfree(ijk)
+      enddo
+
+
+c      scount_buf=scount(me)
+c      scount_buf2=scount(0)
+
+c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
+
 #ifdef MPI
-      call MPI_Gatherv(Fdimless(1),scount(me),
+      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
+      call MPI_Gatherv(Fdimless_buf(1),scount(me),
      & MPI_REAL,Fdimless(1),
      & scount(0),idispl(0),MPI_REAL,Master,
      & MPI_COMM_WORLD, IERROR)
-      call MPI_Gatherv(totfree(1),scount(me),
+      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
+      call MPI_Gatherv(totfree_buf(1),scount(me),
      & MPI_DOUBLE_PRECISION,totfree(1),
      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
      & MPI_COMM_WORLD, IERROR)
-      call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
+      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
+      call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
      & MPI_DOUBLE_PRECISION,entfac(1),
      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
      & MPI_COMM_WORLD, IERROR)
+      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
       if (me.eq.Master) then
+      WRITE (iout,*) "me.eq.Master"
 #endif
 #ifdef DEBUG
         write (iout,*) "The FDIMLESS array before sorting"
@@ -236,13 +260,16 @@ c#endif
 c          write (iout,*) i,fdimless(i)
         enddo
 #endif
+      WRITE (iout,*) "Wchodze do call mysort1"
         call mysort1(ncon,Fdimless,list_conf)
+      WRITE (iout,*) "Wychodze z call mysort1"
 #ifdef DEBUG
         write (iout,*) "The FDIMLESS array after sorting"
         do i=1,ncon
           write (iout,*) i,list_conf(i),fdimless(i)
         enddo
 #endif
+      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
         do i=1,ncon
           totfree(i)=fdimless(i)
         enddo
index 50e1040..f59466c 100644 (file)
@@ -624,9 +624,9 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
         enddo
 #ifndef CLUST
 #ifndef WHAM
-        if (.not.found.and.fg_rank.eq.0) 
-     &      write(iout,'(a15,f12.2,f8.1,2i5)')
-     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+c        if (.not.found.and.fg_rank.eq.0) 
+c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
+c     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
 #endif
 #endif
       enddo
@@ -639,9 +639,9 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
         enddo
 #ifndef CLUST
 #ifndef WHAM
-        if (.not.found.and.fg_rank.eq.0) 
-     &      write(iout,'(a15,f12.2,f8.1,2i5)')
-     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+c        if (.not.found.and.fg_rank.eq.0) 
+c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
+c     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
 #endif
 #endif
       enddo
index cac8518..413817c 100644 (file)
@@ -228,6 +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)
+c     &   +wsccor*fact(1)*gsccor_loc(i)
+c BYLA ROZNICA Z CLUSTER< OSTATNIA LINIA DODANA
       enddo
       endif
       return
@@ -359,11 +361,14 @@ C
       integer icant
       external icant
 cd    print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
+c ROZNICA z cluster
       do i=1,210
         do j=1,2
           eneps_temp(j,i)=0.0d0
         enddo
       enddo
+cROZNICA
+
       evdw=0.0D0
       evdw_t=0.0d0
       do i=iatsc_s,iatsc_e
@@ -397,8 +402,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 cluster
             eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij)
             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)/)')