From 14b55af23c7da34bfdc10f315551c8e01a3a906a Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 19 Apr 2016 08:25:46 +0200 Subject: [PATCH] modyfication for tube --- source/cluster/wham/src-M/probabl.F | 5 +- source/cluster/wham/src-M/read_coords.F | 4 +- source/unres/src_MD-M/MD_A-MTS.F | 1 + source/unres/src_MD-M/energy_p_new_barrier.F | 27 ++- source/unres/src_MD-M/parmread.F | 4 +- source/wham/src-M/enecalc1.F | 4 +- source/wham/src-M/wham_calc1.F | 334 ++++++++++++++++++++++++-- 7 files changed, 338 insertions(+), 41 deletions(-) diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F index 7bc3f1a..67341dd 100644 --- a/source/cluster/wham/src-M/probabl.F +++ b/source/cluster/wham/src-M/probabl.F @@ -296,12 +296,13 @@ c write (iout,*) "qfree",qfree write (iout,*) "ncon", ncon,maxstr_proc do i=1,min0(ncon,maxstr_proc)-1 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree -#ifdef DEBUG +C#ifdef DEBUG + write (iout,*) "tu szukaj ponizej 7" write (iout,*) i,ib,beta_h(ib), & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i), & totfree(list_conf(i)), & -entfac(list_conf(i)),fdimless(i),sumprob -#endif +C#endif if (sumprob.gt.prob_limit) goto 122 c if (sumprob.gt.1.00d0) goto 122 nlist=nlist+1 diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F index 2911e60..55f52e3 100644 --- a/source/cluster/wham/src-M/read_coords.F +++ b/source/cluster/wham/src-M/read_coords.F @@ -393,7 +393,7 @@ c------------------------------------------------------------------------------ chalen=int((nct-nnt+2)/symetr) call int_from_cart1(.false.) do j=nnt+1,nct - if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) + if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.6.0d0) & .and.(itype(j).ne.ntyp1)) then if (j.gt.2) then if (itel(j).ne.0 .and. itel(j-1).ne.0) then @@ -419,7 +419,7 @@ c------------------------------------------------------------------------------ enddo do j=nnt,nct itj=itype(j) - if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0 + if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.5.0d0 & .and. itype(j).ne.ntyp1) then write (iout,*) "Conformation",jjj,jj+1 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 2635a2f..3a45231 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -1653,6 +1653,7 @@ c Removing the velocity of the center of mass #endif call zerograd call etotal(potEcomp) + call enerprint(potEcomp) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 21c0197..f557ebc 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -11536,7 +11536,8 @@ C the vector between center of side_chain and peptide group &pept_group(3),costhet_grad(3),cosphi_grad_long(3), &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3) C the line belowe needs to be changed for FGPROC>1 - do i=1,nres-1 + do i=iatscp_s,iatscp_e +C do i=1,nres-1 if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle ishield_list(i)=0 Cif there two consequtive dummy atoms there is no peptide group between them @@ -11729,8 +11730,8 @@ C lets ommit dummy atoms for now C now calculate distance from center of tube and direction vectors vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize - vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -11779,8 +11780,8 @@ C .or.(iti.eq.10) vectube(1)=mod(vectube(1),boxxsize) if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize vectube(2)=c(2,i+nres) - vectube(2)=mod(vectube(2),boxxsize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(2)=mod(vectube(2),boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -11864,8 +11865,8 @@ C lets ommit dummy atoms for now C now calculate distance from center of tube and direction vectors vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize - vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -11914,8 +11915,8 @@ C in UNRES uncomment the line below as GLY has no side-chain... vectube(1)=mod(vectube(1),boxxsize) if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize vectube(2)=c(2,i+nres) - vectube(2)=mod(vectube(2),boxxsize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize + vectube(2)=mod(vectube(2),boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -11937,10 +11938,10 @@ C lipbufthick is thickenes of lipid buffore ssgradtube=-sscagradlip(fracinbuf)/tubebufthick print *,ssgradtube, sstube,tubetranene(itype(i)) enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) - gg_tube_SC(3,i)=gg_tube_SC(3,i) - &+ssgradtube*tubetranene(itype(i)) - gg_tube(3,i-1)= gg_tube(3,i-1) - &+ssgradtube*tubetranene(itype(i)) +C gg_tube_SC(3,i)=gg_tube_SC(3,i) +C &+ssgradtube*tubetranene(itype(i)) +C gg_tube(3,i-1)= gg_tube(3,i-1) +C &+ssgradtube*tubetranene(itype(i)) C print *,"doing sccale for lower part" elseif (positi.gt.buftubetop) then fracinbuf=1.0d0- diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 5001b6b..8fa2bad 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -1361,7 +1361,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) enddo enddo write(iout,*) "tube param" - read(itube,*) epspeptube,sigmapeptube + read(itube,*) epspeptube,sigmapeptube,tubetranenepep sigmapeptube=sigmapeptube**6 sigeps=dsign(1.0D0,epspeptube) epspeptube=dabs(epspeptube) @@ -1369,7 +1369,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep do i=1,ntyp - read(itube,*) epssctube,sigmasctube + read(itube,*) epssctube,sigmasctube,tubetranene(i) sigmasctube=sigmasctube**6 sigeps=dsign(1.0D0,epssctube) epssctube=dabs(epssctube) diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F index a5d25b3..3a5dbf3 100644 --- a/source/wham/src-M/enecalc1.F +++ b/source/wham/src-M/enecalc1.F @@ -744,7 +744,7 @@ c------------------------------------------------------------------------------ call int_from_cart1(.false.) do j=nnt+1,nct if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. - & (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then + & (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.0d0)) then if (iprint.gt.0) & write (iout,*) "Bad CA-CA bond length",j," ",vbld(j), & " for conformation",ii @@ -769,7 +769,7 @@ c------------------------------------------------------------------------------ do j=nnt,nct itj=itype(j) if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. - & (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then + & (vbld(nres+j)-dsc(iabs(itj))).gt.5.0d0) then if (iprint.gt.0) & write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j), & " for conformation",ii diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index f77a245..f0180ba 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -51,7 +51,8 @@ c parameter (MaxHdim=200000) double precision energia(0:max_ene) #ifdef MPI integer tmax_t,upindE_p - double precision fi_p(MaxR,MaxT_h,Max_Parm) + double precision fi_p(MaxR,MaxT_h,Max_Parm), + & fi_p_min(MaxR,MaxT_h,Max_Parm) double precision sumW_p(0:nGridT,Max_Parm), & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm), & sumQ_p(MaxQ1,0:nGridT,Max_Parm), @@ -63,7 +64,8 @@ c parameter (MaxHdim=200000) & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH, & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h) double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t - double precision potEmin_t,entmin_p,entmax_p + double precision potEmin_t,entmin_p,entmax_p, + & potEmin_t_all(maxT_h,Max_Parm) integer histent_p(0:2000) logical lprint /.true./ #endif @@ -75,11 +77,12 @@ c parameter (MaxHdim=200000) & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT, & weight,econstr double precision fi(MaxR,maxT_h,Max_Parm), + & fi_min(MaxR,maxT_h,Max_Parm), & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom, & hfin(0:MaxHdim,maxT_h),histE(0:maxindE), & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h), - & potEmin,ent, - & hfin_ent(0:MaxHdim),vmax,aux + & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min, + & hfin_ent(0:MaxHdim),vmax,aux,entfac_min double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/, & eplus,eminus,logfac,tanhT,tt @@ -107,6 +110,11 @@ c parameter (MaxHdim=200000) dmin=0.0d0 tmax=0 potEmin=1.0d10 + do i=1,nParmset + do j=1,nT_h(i) + potEmin_all(j,i)=1.0d10 + enddo + enddo rgymin=1.0d10 rmsmin=1.0d10 rgymax=0.0d0 @@ -120,9 +128,9 @@ C#define DEBUG #else do i=1,ntot(islice) #endif - do j=1,nParmSet - if (potE(i,j).le.potEmin) potEmin=potE(i,j) - enddo +C do j=1,nParmSet +C if (potE(i,j).le.potEmin) potEmin=potE(i,j) +C enddo if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i) if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i) if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i) @@ -187,8 +195,8 @@ c & ind_point(i) call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX, & WHAM_COMM,IERROR) tmax=tmax_t - call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION, - & MPI_MIN,WHAM_COMM,IERROR) +C call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION, +C & MPI_MIN,WHAM_COMM,IERROR) call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION, & MPI_MIN,WHAM_COMM,IERROR) call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION, @@ -197,7 +205,8 @@ c & ind_point(i) & MPI_MIN,WHAM_COMM,IERROR) call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION, & MPI_MAX,WHAM_COMM,IERROR) - potEmin=potEmin_t/2 +C potEmin=potEmin_t/2 + rgymin=rgymin_t rgymax=rgymax_t rmsmin=rmsmin_t @@ -392,7 +401,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & *(dd-q0(j,kk,ib,iparm))**2 enddo v(i,kk,ib,iparm)= - & -beta_h(ib,iparm)*(etot-potEmin+Econstr) + & -beta_h(ib,iparm)*(etot+Econstr) #ifdef DEBUG write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm) @@ -408,6 +417,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft ! Compute new free-energy values corresponding to the righ-hand side of the ! equation and their derivatives. write (iout,*) "------------------------fi" + entfac_min=1.0d10 + #ifdef MPI do t=1,scount(me1) #else @@ -433,10 +444,52 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft enddo enddo entfac(t)=-dlog(denom)-vmax + if (entfac(t).lt.entfac_min) entfac_min=entfac(t) #ifdef DEBUG write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t) #endif enddo + + do iparm=1,nParmSet + do iib=1,nT_h(iparm) + do ii=1,nR(iib,iparm) +#ifdef MPI + fi_p_min(ii,iib,iparm)=-1.0d10 + do t=1,scount(me) + aux=v(t,ii,iib,iparm)+entfac(t) + if (aux.gt.fi_p_min(ii,iib,iparm)) + & fi_p_min(ii,iib,iparm)=aux + enddo +#else + do t=1,ntot(islice) + aux=v(t,ii,iib,iparm)+entfac(t) + if (aux.gt.fi_min(ii,iib,iparm)) + & fi_min(ii,iib,iparm)=aux + enddo +#endif + enddo ! ii + enddo ! iib + enddo ! iparm +#ifdef MPI +#ifdef DEBUG + write (iout,*) "fi_min before AllReduce" + do i=1,nParmSet + do j=1,nT_h(i) + write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i)) + enddo + enddo +#endif + call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, + & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR) +#ifdef DEBUG + write (iout,*) "fi_min after AllReduce" + do i=1,nParmSet + do j=1,nT_h(i) + write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i)) + enddo + enddo +#endif +#endif do iparm=1,nParmSet do iib=1,nT_h(iparm) do ii=1,nR(iib,iparm) @@ -444,23 +497,23 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft fi_p(ii,iib,iparm)=0.0d0 do t=1,scount(me) fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) - & +dexp(v(t,ii,iib,iparm)+entfac(t)) + & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm)) #ifdef DEBUG - write (iout,'(4i5,3e15.5)') t,ii,iib,iparm, - & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm) + write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, + & v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), + & fi_p(ii,iib,iparm) #endif enddo #else fi(ii,iib,iparm)=0.0d0 do t=1,ntot(islice) fi(ii,iib,iparm)=fi(ii,iib,iparm) - & +dexp(v(t,ii,iib,iparm)+entfac(t)) + & +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm)) enddo #endif enddo ! ii enddo ! iib enddo ! iparm - #ifdef MPI #ifdef DEBUG write (iout,*) "fi before MPI_Reduce me",me,' master',master @@ -495,7 +548,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft do iparm=1,nParmSet do ib=1,nT_h(iparm) do i=1,nR(ib,iparm) - fi(i,ib,iparm)=-dlog(fi(i,ib,iparm)) + fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm) avefi=avefi+fi(i,ib,iparm) enddo enddo @@ -543,6 +596,198 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft enddo ! iter 20 continue + +#ifdef MPI + do i=1,scount(me1) +#else + do i=1,ntot(islice) +#endif +c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) + do iparm=1,nParmSet +#ifdef DEBUG + write (iout,'(2i5,21f8.2)') i,iparm, + & (enetb(k,i,iparm),k=1,21) +#endif + call restore_parm(iparm) +#ifdef DEBUG + write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, + & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, + & wtor_d,wsccor,wbond +#endif + do ib=1,nT_h(iparm) + if (rescale_mode.eq.1) then + quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) + quotl=1.0d0 + kfacl=1.0d0 + do l=1,5 + quotl1=quotl + quotl=quotl*quot + kfacl=kfacl*kfac + fT(l)=kfacl/(kfacl-1.0d0+quotl) + enddo +#if defined(FUNCTH) + tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3) + ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 +#elif defined(FUNCT) + ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) +#else + ft(6)=1.0d0 +#endif + else if (rescale_mode.eq.2) then + quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) + quotl=1.0d0 + do l=1,5 + quotl=quotl*quot + fT(l)=1.12692801104297249644d0/ + & dlog(dexp(quotl)+dexp(-quotl)) + enddo +#if defined(FUNCTH) + tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3) + ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 +#elif defined(FUNCT) + ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) +#else + ft(6)=1.0d0 +#endif +c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft + else if (rescale_mode.eq.0) then + do l=1,6 + fT(l)=1.0d0 + enddo + else + write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", + & rescale_mode + call flush(iout) + return1 + endif + evdw=enetb(21,i,iparm) + evdw_t=enetb(1,i,iparm) +#ifdef SCP14 + evdw2_14=enetb(17,i,iparm) + evdw2=enetb(2,i,iparm)+evdw2_14 +#else + evdw2=enetb(2,i,iparm) + evdw2_14=0.0d0 +#endif +#ifdef SPLITELE + ees=enetb(3,i,iparm) + evdw1=enetb(16,i,iparm) +#else + ees=enetb(3,i,iparm) + evdw1=0.0d0 +#endif + ecorr=enetb(4,i,iparm) + ecorr5=enetb(5,i,iparm) + ecorr6=enetb(6,i,iparm) + eel_loc=enetb(7,i,iparm) + eello_turn3=enetb(8,i,iparm) + eello_turn4=enetb(9,i,iparm) + eturn6=enetb(10,i,iparm) + ebe=enetb(11,i,iparm) + escloc=enetb(12,i,iparm) + etors=enetb(13,i,iparm) + etors_d=enetb(14,i,iparm) + ehpb=enetb(15,i,iparm) + estr=enetb(18,i,iparm) + esccor=enetb(19,i,iparm) + edihcnstr=enetb(20,i,iparm) +C edihcnstr=0.0d0 + eliptran=enetb(22,i,iparm) +#ifdef SPLITELE + if (shield_mode.gt.0) then + etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 + & +ft(1)*welec*ees + & +ft(1)*wvdwpp*evdw1 + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +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+wliptran*eliptran + else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees + & +wvdwpp*evdw1 + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +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+wliptran*eliptran + endif +#else + if (shield_mode.gt.0) then + etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 + & +ft(1)*welec*(ees+evdw1) + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +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+wliptran*eliptran + else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 + & +ft(1)*welec*(ees+evdw1) + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +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+wliptran*eliptran + endif + +#endif + etot=etot-entfac(i)/beta_h(ib,iparm) + if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot + + enddo ! ib + enddo ! iparm + enddo ! i + + +#ifdef DEBUG + write (iout,*) "The potEmin array before reduction" + do i=1,nParmSet + write (iout,*) "Parameter set",i + do j=1,nT_h(i) + write (iout,*) j,PotEmin_all(j,i) + enddo + enddo + write (iout,*) "potEmin_min",potEmin_min +#endif +#ifdef MPI +C Determine the minimum energes for all parameter sets and temperatures + call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), + & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR) + do i=1,nParmSet + do j=1,nT_h(i) + potEmin_all(j,i)=potEmin_t_all(j,i) + enddo + enddo +#endif + potEmin_min=potEmin_all(1,1) + do i=1,nParmSet + do j=1,nT_h(i) + if (potEmin_all(j,i).lt.potEmin_min) + & potEmin_min=potEmin_all(j,i) + enddo + enddo +#ifdef DEBUG + write (iout,*) "The potEmin array" + do i=1,nParmSet + write (iout,*) "Parameter set",i + do j=1,nT_h(i) + write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), + & PotEmin_all(j,i) + enddo + enddo + write (iout,*) "potEmin_min",potEmin_min +#endif + + ! Now, put together the histograms from all simulations, in order to get the ! unbiased total histogram. #ifdef MPI @@ -717,7 +962,9 @@ c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18) estr=enetb(18,t,iparm) esccor=enetb(19,t,iparm) edihcnstr=enetb(20,t,iparm) - edihcnstr=0.0d0 +C edihcnstr=0.0d0 + eliptran=enetb(22,i,iparm) + do k=0,nGridT betaT=startGridT+k*delta_T temper=betaT @@ -798,6 +1045,36 @@ c ft=2*T0/(T0+betaT) c write (iout,*) "ftprim",ftprim c write (iout,*) "ftbis",ftbis betaT=1.0d0/(1.987D-3*betaT) + if (betaT.ge.beta_h(1,iparm)) then + potEmin=potEmin_all(1,iparm)+ + & (potEmin_all(1,iparm)-potEmin_all(2,iparm))/ + & (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))* + & (1.0/betaT-1.0/beta_h(1,iparm)) +#ifdef DEBUG + write(iout,*) "first",temper,potEmin +#endif + else if (betaT.le.beta_h(nT_h(iparm),iparm)) then + potEmin=potEmin_all(nT_h(iparm),iparm)+ + &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/ + &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))* + &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm)) +#ifdef DEBUG + write (iout,*) "last",temper,potEmin +#endif + else + do l=1,nT_h(iparm)-1 + if (betaT.le.beta_h(l,iparm) .and. + & betaT.gt.beta_h(l+1,iparm)) then + potEmin=potEmin_all(l,iparm) +#ifdef DEBUG + write (iout,*) "l",l, + & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)), + & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin +#endif + exit + endif + enddo + endif #ifdef SPLITELE if (shield_mode.gt.0) then etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 @@ -882,7 +1159,7 @@ C & +ftprim(6)*evdw_t & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ - & ftprim(1)*wsccor*esccor + & ftbis(1)*wsccor*esccor else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 & +ft(1)*welec*(ees+evdw1) @@ -905,7 +1182,7 @@ C & +ftprim(6)*evdw_t & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ - & ftprim(1)*wsccor*esccor + & ftbis(1)*wsccor*esccor endif @@ -948,6 +1225,7 @@ C & +ftprim(6)*evdw_t endif #ifdef MPI do ib=1,nT_h(iparm) + potEmin=potEmin_all(ib,iparm) expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) hfin_p(ind,ib)=hfin_p(ind,ib)+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) @@ -960,6 +1238,7 @@ C & +ftprim(6)*evdw_t enddo #else do ib=1,nT_h(iparm) + potEmin=potEmin_all(ib,iparm) expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) hfin(ind,ib)=hfin(ind,ib)+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) @@ -1178,6 +1457,21 @@ C & +ftprim(6)*evdw_t write (iout,'(a,i3)') "Parameter set",iparm endif do i=0,NGridT + betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T)) + if (betaT.ge.beta_h(1,iparm)) then + potEmin=potEmin_all(1,iparm) + else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then + potEmin=potEmin_all(nT_h(iparm),iparm) + else + do l=1,nT_h(iparm)-1 + if (betaT.le.beta_h(l,iparm) .and. + & betaT.gt.beta_h(l+1,iparm)) then + potEmin=potEmin_all(l,iparm) + exit + endif + enddo + endif + sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm) sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ & sumW(i,iparm) -- 1.7.9.5