Fixed main module of WHAM

diff --git a/source/wham/src/Makefile_MPICH_ifort b/source/wham/src/Makefile_MPICH_ifort
index 6e2ba17..711b6c7 100644 (file)
--- a/source/wham/src/Makefile_MPICH_ifort
+++ b/source/wham/src/Makefile_MPICH_ifort
@@ -71,7 +71,7 @@ E0LL2Y: ${objects} ${objects_compar} xdrf/libxdrf.a
        ./compinfo
        ${FC} -c ${FFLAGS} cinfo.f
        $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \
-       ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH-restr-DFA_E0LL2Y.exe
+       ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH-restr-DFA_E0LL2Y-test.exe
 
 xdrf/libxdrf.a:
        cd xdrf && make

diff --git a/source/wham/src/cxread.F b/source/wham/src/cxread.F
index a662f7a..7c7ae50 100644 (file)
--- a/source/wham/src/cxread.F
+++ b/source/wham/src/cxread.F
@@ -5,6 +5,7 @@
       include 'DIMENSIONS.FREE'
       integer MaxTraj
       parameter (MaxTraj=2050)
+      include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
@@ -81,7 +82,8 @@ c      print *,"bumbum"
        endif
       enddo
       call xdrfint_(ixdrf, nprop, iret)
-      if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep) 
+      if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
+     &  .or. hamil_rep) 
      &  call xdrfint(ixdrf, iset, iret)
       do i=1,nprop
         call xdrffloat_(ixdrf, rprop(i), iret)
@@ -124,7 +126,8 @@ c      write (iout,*) "nprop",nprop
         nprop_prev=nprop
       endif
       call flush(iout)
-      if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep) 
+      if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
+     &  .or. hamil_rep) 
      &  call xdrfint(ixdrf, iset, iret)
       do i=1,nprop
         call xdrffloat(ixdrf, rprop(i), iret)
@@ -254,7 +257,7 @@ c            write (iout,*) "Opening file ",
 c     &          bprotfile_temp(:ilen(bprotfile_temp))
             islice1=islice
           endif
-          if (umbrella(iparm)) then
+          if (umbrella(iparm) .or. homol_nset.gt.1) then
             write(ientout,rec=ll(islice))
      &        ((xoord(l,k),l=1,3),k=1,nres),
      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
@@ -277,6 +280,8 @@ c     &          bprotfile_temp(:ilen(bprotfile_temp))
      &        iR,iib,iparm
           endif
 #ifdef DEBUG
+          write (iout,*) " constr_homology",constr_homology,
+     &      " ll",ll(islice)," iset",iset
           call int_from_cart1(.false.)
           write (iout,*) "Writing conformation, record",ll(islice)
           write (iout,*) "Cartesian coordinates"

diff --git a/source/wham/src/enecalc1.F b/source/wham/src/enecalc1.F
index dcf6a97..bdda141 100644 (file)
--- a/source/wham/src/enecalc1.F
+++ b/source/wham/src/enecalc1.F
@@ -77,6 +77,15 @@
      &    ((csingle(l,k+nres),l=1,3),k=nnt,nct),
      &    nss,(ihpb(k),jhpb(k),k=1,nss),
      &    eini,efree,rmsdev,(q(j,iii+1),j=1,nQ),iR,ib,ipar
+#ifdef DEBUG
+        write (iout,*) "homol_nset",homol_nset,
+     &    "  i",i," iR",iR," ib",ib," iset",iset
+#endif
+        if (homol_nset.gt.1) iset=iR
+#ifdef DEBUG
+        write (iout,*) "homol_nset",homol_nset,
+     &    "  i",i," iR",iR," ib",ib," iset",iset
+#endif
 cc       write(iout,*), 'NAWEJ',i,eini
          if (indpdb.gt.0) then
            do k=1,nres
@@ -156,6 +165,8 @@ c     &   " kfac",kfac,"quot",quot," fT",fT
      &      wtor_d,wsccor,wbond
 #endif
         call etotal(energia(0),fT)
+        if (constr_homology) energia(0)=energia(0)+
+     &   waga_homology(iset)*energia(22)
 #ifdef DEBUG
         write (iout,*) "Conformation",i
         call enerprint(energia(0),fT)

diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F
index 3ca5082..52bc9a6 100644 (file)
--- a/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@ -133,7 +133,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
+     & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr
      & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
      & +wdfa_beta*edfabet
 #else
@@ -144,7 +144,7 @@ c      write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t
      & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4
      & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
      & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
-     & +wbond*estr+wsccor*fact(1)*esccor+ehomology_constr
+     & +wbond*estr+wsccor*fact(1)*esccor!+ehomology_constr
      & +wdfa_dist*edfadis+wdfa_tor*edfator+wdfa_nei*edfanei
      & +wdfa_beta*edfabet
 #endif
@@ -3706,15 +3706,19 @@ c
 c
 c          For Gaussian-type Urestr
 c
-        ehomology_constr=(waga_dist*odleg+waga_angle*kat+
-     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c        ehomology_constr=(waga_dist*odleg+waga_angle*kat+
+c     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+        ehomology_constr=waga_dist*odleg+waga_angle*kat+
+     &              waga_theta*Eval+waga_d*Erot
 c     write (iout,*) "ehomology_constr=",ehomology_constr
       else
 c
 c          For Lorentzian-type Urestr
 c  
-        ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
-     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+c        ehomology_constr=(-waga_dist*odleg+waga_angle*kat+
+c     &              waga_theta*Eval+waga_d*Erot)*waga_homology(iset)
+        ehomology_constr=-waga_dist*odleg+waga_angle*kat+
+     &              waga_theta*Eval+waga_d*Erot
 c     write (iout,*) "ehomology_constr=",ehomology_constr
       endif
 #ifdef DEBUG

diff --git a/source/wham/src/readrtns.F b/source/wham/src/readrtns.F
index ad038e2..19d4b20 100644 (file)
--- a/source/wham/src/readrtns.F
+++ b/source/wham/src/readrtns.F
@@ -147,6 +147,11 @@ C
       endif
       replica(iparm)=index(controlcard,"REPLICA").gt.0
       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
+      if (umbrella(iparm) .and. homol_nset.gt.1) then
+        umbrella(iparm) = .false.
+        write (iout,*) 
+     & "Replica in homology restraints weights UMBRELLA ignored,",iparm
+      endif
       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
      &  replica(iparm)," umbrella ",umbrella(iparm),

diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F
index d9c3de2..986f6a3 100644 (file)
--- a/source/wham/src/wham_calc1.F
+++ b/source/wham/src/wham_calc1.F
@@ -328,7 +328,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
+     &      +wbond*estr+wdfa_dist*edfadis
      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 #else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
@@ -339,7 +339,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
+     &      +wbond*estr+wdfa_dist*edfadis
      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 #endif
 #ifdef DEBUG
@@ -356,6 +356,21 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             call enerprint(energia(0),fT)
             endif
 #endif
+            if (homol_nset.le.1) then
+
+            do kk=1,nR(ib,iparm)
+              Econstr=waga_homology(kk)*ehomology_constr
+              v(i,kk,ib,iparm)=
+     &          -beta_h(ib,iparm)*(etot+Econstr)
+#ifdef DEBUG
+              write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
+     &         etot,v(i,kk,ib,iparm)
+#endif
+            enddo ! kk
+
+            else
+
+            etot=etot+ehomology_constr
             do kk=1,nR(ib,iparm)
               Econstr=0.0d0
               do j=1,nQ
@@ -370,6 +385,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &         etot,v(i,kk,ib,iparm)
 #endif
             enddo ! kk
+
+            endif
           enddo   ! ib
         enddo     ! iparm
       enddo       ! i