changes need for chamm-gui
[unres4.git] / source / wham / enecalc.F90
index 1568a9a..6ba59a0 100644 (file)
       character(len=3) :: liczba
 !el      real(kind=8) :: qwolynes
 !el      external qwolynes
-      integer :: errmsg_count,maxerrmsg_count=100
+      integer :: errmsg_count,maxerrmsg_count=100000
 !el      real(kind=8) :: rmsnat,gyrate
 !el      external rmsnat,gyrate
-      real(kind=8) :: tole=1.0d-1
+      real(kind=8) :: tole=0.0d0
       integer i,itj,ii,iii,j,k,l,licz
       integer ir,ib,ipar,iparm
       integer iscor,islice
@@ -405,7 +405,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot
             endif
           endif
           potE(iii+1,iparm)=energia(0)
-          do k=1,49
+          do k=1,max_ene
             enetb(k,iii+1,iparm)=energia(k)
           enddo
 !           write (iout,'(2i5,21f8.2)') "debug",k,iii+1,(enetb(k,iii+1,iparm),k=1,21)
@@ -1231,7 +1231,9 @@ write(iout,*)"end of store_parm"
             eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, &
             ecation_prot, ecationcation,evdwpp,eespp,evdwpsb,eelpsb,evdwsb, &
             eelsb,estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
-            ecorr_nucl,ecorr3_nucl,escbase, epepbase,escpho, epeppho
+            ecorr_nucl,ecorr3_nucl,escbase, epepbase,escpho, epeppho,ecation_nucl,&
+            elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang, &
+            eliptran,ehomology_constr
 
       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
@@ -1246,7 +1248,7 @@ write(iout,*)"end of store_parm"
       integer :: iperm(MaxStr)
       integer :: islice
       integer,dimension(0:nprocs) :: scount_
+      write(iout,*) "Begin make ensemble" 
 #ifdef MPI
       if (me.eq.Master) then
 #endif
@@ -1271,14 +1273,15 @@ write(iout,*)"end of store_parm"
         form="unformatted",access="direct",recl=lenrec1)
 #ifdef MPI
       endif
-#endif
+#endif 
+      write(iout,*) "iparmprint",iparmprint,iparm
       do iparm=1,iparm
         if (iparm.ne.iparmprint) exit
         call restore_parm(iparm)
         do ib=1,nT_h(iparm)
-#ifdef DEBUG
+!#ifdef DEBUG
           write (iout,*) "iparm",iparm," ib",ib
-#endif
+!#endif
           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
 !          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 !          quotl=1.0d0
@@ -1396,6 +1399,25 @@ write(iout,*)"end of store_parm"
             etors_d_nucl =      enetb(36,i,iparm)
             ecorr_nucl =      enetb(37,i,iparm)
             ecorr3_nucl =      enetb(38,i,iparm)
+            epeppho=   enetb(49,i,iparm)
+            escpho=    enetb(48,i,iparm)
+            epepbase=  enetb(47,i,iparm)
+            escbase=   enetb(46,i,iparm)
+            ecation_nucl=enetb(50,i,iparm)
+            elipbond=enetb(52,i,iparm)
+            elipang=enetb(53,i,iparm)
+            eliplj=enetb(54,i,iparm)
+            elipelec=enetb(55,i,iparm)
+            ecat_prottran=enetb(56,i,iparm)
+            ecation_protang=enetb(57,i,iparm)
+            eliptran=enetb(22,i,iparm)
+            ehomology_constr=enetb(51,i,iparm)
+          if (homol_nset.gt.1) &
+            ehomology_constr=waga_homology(homol_nset)*ehomology_constr
+!      wscbase=ww(46)
+!      wpepbase=ww(47)
+!      wscpho=ww(48)
+!      wpeppho=ww(49)
 
 !#ifdef SPLITELE
 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
@@ -1423,7 +1445,7 @@ write(iout,*)"end of store_parm"
             etot=wsc*evdw+wscp*evdw2+welec*ees &
             +wvdwpp*evdw1 &
             +wang*ebe+wtor*etors+wscloc*escloc &
-            +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+            +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
             +wcorr6*ecorr6+wturn4*eello_turn4 &
             +wturn3*eello_turn3 &
             +wturn6*eello_turn6+wel_loc*eel_loc &
@@ -1433,14 +1455,18 @@ write(iout,*)"end of store_parm"
             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
-            +wscbase*escbase&
-            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+            +wscbase*escbase+ehomology_constr&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+            +wcatnucl*ecation_nucl&
+            +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang+&
+             wliptran*eliptran
+
 
 #else
             etot=wsc*evdw+wscp*evdw2 &
             +welec*(ees+evdw1) &
             +wang*ebe+wtor*etors+wscloc*escloc &
-            +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+            +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
             +wcorr6*ecorr6+wturn4*eello_turn4 &
             +wturn3*eello_turn3 &
             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
@@ -1450,8 +1476,13 @@ write(iout,*)"end of store_parm"
             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
-            +wscbase*escbase&
-            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+            +wscbase*escbase+ehomology_constr&
+            +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+            +wcatnucl*ecation_nucl&
+            +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+            +wliptran*eliptran
+
+
 
 #endif
 
@@ -1459,24 +1490,33 @@ write(iout,*)"end of store_parm"
             Fdimless(i)= &
               beta_h(ib,iparm)*etot-entfac(i)
             potE(i,iparm)=etot
-#ifdef DEBUG
+!#ifdef DEBUG
             write (iout,*) i,indstart(me)+i-1,ib,&
              1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),&
              -entfac(i),Fdimless(i)
-#endif
+!#endif
 #else
             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
             potE(i,iparm)=etot
 #endif
           enddo   ! i
 #ifdef MPI
+          scount_(:)=0
           do i=1,scount(me1)
             Fdimless_(i)=Fdimless(i)
           enddo
+          write(iout,*) "before gather Fdimless"
+          write(iout,*) scount(me),scount_(0),idispl(0)
+          write (iout,*) "added update of scount_"
+         call MPI_AllGather( scount(me), 1, MPI_INTEGER, scount_(0), 1, &
+        MPI_INTEGER, WHAM_COMM, IERROR)
+
+
           call MPI_Gatherv(Fdimless_(1),scount(me),&
            MPI_REAL,Fdimless(1),&
            scount_(0),idispl(0),MPI_REAL,Master,&
            WHAM_COMM, IERROR)
+          write(iout,*) "after gather Fdimless"
 #ifdef DEBUG
           call MPI_Gatherv(potE(1,iparm),scount_(me),&
            MPI_DOUBLE_PRECISION,potE(1,iparm),&
@@ -1488,23 +1528,23 @@ write(iout,*)"end of store_parm"
            WHAM_COMM, IERROR)
 #endif
           if (me.eq.Master) then
-#ifdef DEBUG
+!#ifdef DEBUG
           write (iout,*) "The FDIMLESS array before sorting"
           do i=1,ntot(islice)
             write (iout,*) i,fdimless(i)
           enddo
-#endif
+!#endif
 #endif
           do i=1,ntot(islice)
             iperm(i)=i
           enddo
           call mysort1(ntot(islice),Fdimless,iperm)
-#ifdef DEBUG
+!#ifdef DEBUG
           write (iout,*) "The FDIMLESS array after sorting"
           do i=1,ntot(islice)
             write (iout,*) i,iperm(i),fdimless(i)
           enddo
-#endif
+!#endif
           qfree=0.0d0
           do i=1,ntot(islice)
             qfree=qfree+exp(-fdimless(i)+fdimless(1))
@@ -1514,12 +1554,12 @@ write(iout,*)"end of store_parm"
           sumprob=0.0
           do i=1,min0(ntot(islice),ensembles) 
             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
-#ifdef DEBUG
+!#ifdef DEBUG
             write (iout,*) i,ib,beta_h(ib,iparm),&
              1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),&
              potE(iperm(i),iparm),&
              -entfac(iperm(i)),fdimless(i),sumprob
-#endif
+!#endif
             if (sumprob.gt.0.99d0) goto 122
             nlist=nlist+1
           enddo  
@@ -1549,10 +1589,10 @@ write(iout,*)"end of store_parm"
             ik=ii-indstart(iproc)+1
             if (iproc.ne.Master) then
               if (me.eq.iproc) then
-#ifdef DEBUG
+!#ifdef DEBUG
                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,&
                  " energy",potE(ik,iparm)
-#endif
+!#endif
                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,&
                   Master,i,WHAM_COMM,IERROR)
               else if (me.eq.Master) then
@@ -1568,6 +1608,7 @@ write(iout,*)"end of store_parm"
             enepot(i)=potE(iperm(i),iparm)
           enddo
 #endif
+          write(iout,*) "DEBUG",me
 #ifdef MPI
           if (me.eq.Master) then
 #endif