1 subroutine make_ensembles(islice,*)
2 ! construct the conformational ensembles at REMD temperatures
5 include "DIMENSIONS.ZSCOPT"
6 include "DIMENSIONS.FREE"
10 integer ierror,errcode,status(MPI_STATUS_SIZE)
12 include "COMMON.IOUNITS"
13 include "COMMON.CONTROL"
15 include "COMMON.ENERGIES"
16 include "COMMON.FFIELD"
17 include "COMMON.INTERACT"
18 include "COMMON.SBRIDGE"
19 include "COMMON.CHAIN"
20 include "COMMON.PROTFILES"
22 real*4 csingle(3,maxres2)
23 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
24 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
25 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
26 & escloc,ehomology_constr,esaxs,
27 & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
28 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
29 integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
30 double precision qfree,sumprob,eini,efree,rmsdev
32 character*2 licz1,licz2
33 character*3 licz3,licz4
37 real*4 Fdimless(MaxStr),Fdimless_(MaxStr)
38 double precision enepot(MaxStr)
43 if (me.eq.Master) then
45 write (licz2,'(bz,i2.2)') islice
47 if (.not.separate_parset) then
48 bxname = prefix(:ilen(prefix))//".bx"
50 write (licz3,'(bz,i3.3)') myparm
51 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
54 if (.not.separate_parset) then
55 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
57 write (licz3,'(bz,i3.3)') myparm
58 bxname = prefix(:ilen(prefix))//"par_"//licz3//
59 & "_slice_"//licz2//".bx"
62 open (ientout,file=bxname,status="unknown",
63 & form="unformatted",access="direct",recl=lenrec1)
68 if (iparm.ne.iparmprint) exit
69 call restore_parm(iparm)
72 write (iout,*) "iparm",iparm," ib",ib
74 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
75 c quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
82 c fT(l)=kfacl/(kfacl-1.0d0+quotl)
84 if (rescale_mode.eq.1) then
85 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
87 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
88 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
100 fT(l)=kfacl/(kfacl-1.0d0+quotl)
102 else if (rescale_mode.eq.2) then
103 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
105 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
106 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
115 fT(l)=1.12692801104297249644d0/
116 & dlog(dexp(quotl)+dexp(-quotl))
118 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
119 else if (rescale_mode.eq.0) then
125 & "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
134 evdw=enetb(1,i,iparm)
135 evdw_t=enetb(21,i,iparm)
137 evdw2_14=enetb(17,i,iparm)
138 evdw2=enetb(2,i,iparm)+evdw2_14
140 evdw2=enetb(2,i,iparm)
145 evdw1=enetb(16,i,iparm)
150 ecorr=enetb(4,i,iparm)
151 ecorr5=enetb(5,i,iparm)
152 ecorr6=enetb(6,i,iparm)
153 eel_loc=enetb(7,i,iparm)
154 eello_turn3=enetb(8,i,iparm)
155 eello_turn4=enetb(9,i,iparm)
156 eturn6=enetb(10,i,iparm)
157 ebe=enetb(11,i,iparm)
158 escloc=enetb(12,i,iparm)
159 etors=enetb(13,i,iparm)
160 etors_d=enetb(14,i,iparm)
161 ehpb=enetb(15,i,iparm)
162 estr=enetb(18,i,iparm)
163 esccor=enetb(19,i,iparm)
164 edihcnstr=enetb(20,i,iparm)
165 ehomology_constr=enetb(22,i,iparm)
166 esaxs=enetb(26,i,iparm)
168 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
170 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
171 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
172 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
173 & +ft(2)*wturn3*eello_turn3
174 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
175 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
176 & +wbond*estr+ehomology_constr+wsaxs*esaxs
178 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
179 & +ft(1)*welec*(ees+evdw1)
180 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
181 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
182 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
183 & +ft(2)*wturn3*eello_turn3
184 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
185 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
186 & +wbond*estr+ehomology_constr+wsaxs*esaxs
190 & beta_h(ib,iparm)*etot-entfac(i)
193 write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
194 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
195 & -entfac(i),Fdimless(i)
198 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
203 call MPI_Gatherv(Fdimless_(1),scount(me),
204 & MPI_REAL,Fdimless(1),
205 & scount(0),idispl(0),MPI_REAL,Master,
208 call MPI_Gatherv(potE(1,iparm),scount(me),
209 & MPI_DOUBLE_PRECISION,potE(1,iparm),
210 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
212 call MPI_Gatherv(entfac(1),scount(me),
213 & MPI_DOUBLE_PRECISION,entfac(1),
214 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
217 if (me.eq.Master) then
219 write (iout,*) "The FDIMLESS array before sorting"
221 write (iout,*) i,fdimless(i)
228 call mysort1(ntot(islice),Fdimless,iperm)
230 write (iout,*) "The FDIMLESS array after sorting"
232 write (iout,*) i,iperm(i),fdimless(i)
237 qfree=qfree+exp(-fdimless(i)+fdimless(1))
239 c write (iout,*) "qfree",qfree
242 do i=1,min0(ntot(islice),ensembles)
243 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
245 write (iout,*) i,ib,beta_h(ib,iparm),
246 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
247 & potE(iperm(i),iparm),
248 & -entfac(iperm(i)),fdimless(i),sumprob
250 if (sumprob.gt.0.99d0) goto 122
256 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
258 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
263 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
266 if (iproc.ge.nprocs) then
267 write (iout,*) "Fatal error: processor out of range",iproc
272 close (ientout,status="delete")
276 ik=ii-indstart(iproc)+1
277 if (iproc.ne.Master) then
278 if (me.eq.iproc) then
280 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
281 & " energy",potE(ik,iparm)
283 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
284 & Master,i,WHAM_COMM,IERROR)
285 else if (me.eq.Master) then
286 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
287 & WHAM_COMM,STATUS,IERROR)
289 else if (me.eq.Master) then
290 enepot(i)=potE(ik,iparm)
295 enepot(i)=potE(iperm(i),iparm)
299 if (me.eq.Master) then
301 write(licz3,'(bz,i3.3)') iparm
302 write(licz2,'(bz,i2.2)') islice
303 if (temper.lt.100.0d0) then
304 write(ctemper,'(f3.0)') temper
305 else if (temper.lt.1000.0) then
306 write (ctemper,'(f4.0)') temper
308 write (ctemper,'(f5.0)') temper
310 if (nparmset.eq.1) then
311 if (separate_parset) then
312 write(licz4,'(bz,i3.3)') myparm
313 pdbname=prefix(:ilen(prefix))//"_par"//licz4
315 pdbname=prefix(:ilen(prefix))
318 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
320 if (nslice.eq.1) then
321 pdbname=pdbname(:ilen(pdbname))//"_T_"//
322 & ctemper(:ilen(ctemper))//"pdb"
324 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
325 & ctemper(:ilen(ctemper))//"pdb"
327 open(ipdb,file=pdbname)
329 read (ientout,rec=iperm(i))
330 & ((csingle(l,k),l=1,3),k=1,nres),
331 & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
332 & nss,(ihpb(k),jhpb(k),k=1,nss),
333 & eini,efree,rmsdev,iscor
340 call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
350 close(ientout,status="delete")
354 !--------------------------------------------------
355 subroutine mysort1(n, x, ipermut)
357 integer i,j,imax,ipm,n
365 if (x(j).lt.xtemp) then
373 ipermut(imax)=ipermut(i)