From 9e7d8385e56076981e143dc87e89d3e4b4a0abe4 Mon Sep 17 00:00:00 2001 From: Adam Liwo Date: Thu, 28 May 2015 07:07:13 +0200 Subject: [PATCH] Fixed main module of WHAM --- source/wham/src/Makefile_MPICH_ifort | 2 +- source/wham/src/cxread.F | 11 ++++++++--- source/wham/src/enecalc1.F | 11 +++++++++++ source/wham/src/energy_p_new.F | 16 ++++++++++------ source/wham/src/readrtns.F | 5 +++++ source/wham/src/wham_calc1.F | 21 +++++++++++++++++++-- 6 files changed, 54 insertions(+), 12 deletions(-) diff --git a/source/wham/src/Makefile_MPICH_ifort b/source/wham/src/Makefile_MPICH_ifort index 6e2ba17..711b6c7 100644 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 -- 1.7.9.5