X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fenecalc.F90;h=d0bd838d476a32a808cd7436bc6a54fc9261289b;hb=bc23440fbe68672d430f71f22f46b11265f003db;hp=b8fd48ff382660148d47e2190fedaeda4e0d55cf;hpb=df2469d9ac903d93889867f4e50e9bf6c428c1c6;p=unres4.git diff --git a/source/wham/enecalc.F90 b/source/wham/enecalc.F90 index b8fd48f..d0bd838 100644 --- a/source/wham/enecalc.F90 +++ b/source/wham/enecalc.F90 @@ -2,6 +2,7 @@ !----------------------------------------------------------------------------- use io_units use wham_data + use control_data, only: tor_mode ! use geometry_data, only:nres,boxxsize,boxysize,boxzsize use energy_data @@ -404,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,50 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) @@ -555,7 +556,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then if (iprint.gt.0) & write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),& - " for conformation",ii + " for conformation",ii,mnum if (iprint.gt.1) then write (iout,*) "The Cartesian geometry is:" write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) @@ -675,6 +676,8 @@ write(iout,*)"enecalc_ i ntot",i,ntot ww_all(16,iparm)=wvdwpp ww_all(17,iparm)=wbond ww_all(19,iparm)=wsccor + ww_all(42,iparm)=wcatprot + ww_all(41,iparm)=wcatcat ! Store bond parameters vbldp0_all(iparm)=vbldp0 akp_all(iparm)=akp @@ -709,6 +712,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot sigc0_all(i,iparm)=sigc0(i) enddo #else + IF (tor_mode.eq.0) THEN nthetyp_all(iparm)=nthetyp ntheterm_all(iparm)=ntheterm ntheterm2_all(iparm)=ntheterm2 @@ -760,6 +764,9 @@ write(iout,*)"enecalc_ i ntot",i,ntot enddo enddo enddo + ELSE + write(iout,*) "Need storing parameters" + ENDIF #endif #ifdef CRYST_SC ! Store the sidechain rotamer parameters @@ -787,6 +794,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot enddo enddo #endif + IF (tor_mode.eq.0) THEN ! Store the torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 @@ -829,7 +837,7 @@ write(iout,*)"enecalc_ i ntot",i,ntot enddo enddo ! Store parameters of the cumulants - do i=-nloctyp,nloctyp + do i=1,nloctyp do j=1,2 b1_all(j,i,iparm)=b1(j,i) b1tilde_all(j,i,iparm)=b1tilde(j,i) @@ -845,6 +853,9 @@ write(iout,*)"enecalc_ i ntot",i,ntot enddo enddo enddo + ELSE + write(iout,*) "NEED storing parameters" + ENDIF ! Store the parameters of electrostatic interactions do i=1,2 do j=1,2 @@ -947,6 +958,8 @@ write(iout,*)"end of store_parm" wvdwpp=ww_all(16,iparm) wbond=ww_all(17,iparm) wsccor=ww_all(19,iparm) + wcatprot=ww_all(42,iparm) + wcatcat=ww_all(41,iparm) ! Restore bond parameters vbldp0=vbldp0_all(iparm) akp=akp_all(iparm) @@ -981,6 +994,7 @@ write(iout,*)"end of store_parm" sigc0(i)=sigc0_all(i,iparm) enddo #else + IF (tor_mode.eq.0) THEN nthetyp=nthetyp_all(iparm) ntheterm=ntheterm_all(iparm) ntheterm2=ntheterm2_all(iparm) @@ -1032,6 +1046,9 @@ write(iout,*)"end of store_parm" enddo enddo enddo + ELSE + write (iout,*) "Need storing parameters" + ENDIF #endif ! Restore the sidechain rotamer parameters #ifdef CRYST_SC @@ -1058,6 +1075,7 @@ write(iout,*)"end of store_parm" enddo enddo #endif + IF (tor_mode.eq.0) THEN ! Restore the torsional parameters do iblock=1,2 do i=-ntortyp+1,ntortyp-1 @@ -1099,8 +1117,9 @@ write(iout,*)"end of store_parm" enddo enddo enddo + ! Restore parameters of the cumulants - do i=-nloctyp,nloctyp + do i=1,nloctyp do j=1,2 b1(j,i)=b1_all(j,i,iparm) b1tilde(j,i)=b1tilde_all(j,i,iparm) @@ -1116,6 +1135,9 @@ write(iout,*)"end of store_parm" enddo enddo enddo + ELSE + write(iout,*) "need storing parameters" + ENDIF ! Restore the parameters of electrostatic interactions do i=1,2 do j=1,2 @@ -1207,7 +1229,10 @@ write(iout,*)"end of store_parm" real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,& eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, & - ecation_prot, ecationcation + 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,ecation_nucl + integer :: i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist real(kind=8) :: qfree,sumprob,eini,efree,rmsdev character(len=80) :: bxname @@ -1221,7 +1246,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 @@ -1246,14 +1271,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 @@ -1358,6 +1384,28 @@ write(iout,*)"end of store_parm" edihcnstr=enetb(19,i,iparm) ecationcation=enetb(42,i,iparm) ecation_prot=enetb(41,i,iparm) + evdwpp = enetb(26,i,iparm) + eespp = enetb(27,i,iparm) + evdwpsb = enetb(28,i,iparm) + eelpsb = enetb(29,i,iparm) + evdwsb = enetb(30,i,iparm) + eelsb = enetb(31,i,iparm) + estr_nucl = enetb(32,i,iparm) + ebe_nucl = enetb(33,i,iparm) + esbloc = enetb(34,i,iparm) + etors_nucl = enetb(35,i,iparm) + 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) +! 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 & @@ -1385,46 +1433,71 @@ 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 & +edihcnstr+wtor_d*etors_d+wsccor*esccor & - +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation& + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& + +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& + +wcatnucl*ecation_nucl + #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 & +wtor_d*etors_d+wsccor*esccor & - +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation + +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation& + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl& + +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& + +wcatnucl*ecation_nucl + #endif #ifdef MPI 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),& @@ -1436,23 +1509,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)) @@ -1462,12 +1535,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 @@ -1497,10 +1570,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 @@ -1516,6 +1589,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