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,eliptran,etube,
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)
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 eliptran=enetb(22,i,iparm)
166 etube=enetb(25,i,iparm)
169 if (shield_mode.gt.0) then
170 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
172 & +ft(1)*wvdwpp*evdw1
173 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
174 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
175 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
176 & +ft(2)*wturn3*eello_turn3
177 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
178 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
179 & +wbond*estr+wliptran*eliptran+wtube*etube
181 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
183 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
184 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
185 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
186 & +ft(2)*wturn3*eello_turn3
187 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
188 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
189 & +wbond*estr+wliptran*eliptran+wtube*etube
192 if (shield_mode.gt.0) then
193 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
194 & +ft(1)*welec*(ees+evdw1)
195 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
196 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
197 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
198 & +ft(2)*wturn3*eello_turn3
199 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
200 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
201 & +wbond*estr+wliptran*eliptran+wtube*etube
203 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
204 & +ft(1)*welec*(ees+evdw1)
205 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
206 & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
207 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
208 & +ft(2)*wturn3*eello_turn3
209 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
210 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
211 & +wbond*estr+wliptran*eliptran+wtube*etube
217 & beta_h(ib,iparm)*etot-entfac(i)
220 write (iout,*) i,indstart(me)+i-1,ib,
221 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
222 & -entfac(i),Fdimless(i)
225 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
230 call MPI_Gatherv(Fdimless(1),scount(me),
231 & MPI_REAL,Fdimless(1),
232 & scount(0),idispl(0),MPI_REAL,Master,
235 call MPI_Gatherv(potE(1,iparm),scount(me),
236 & MPI_DOUBLE_PRECISION,potE(1,iparm),
237 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
239 call MPI_Gatherv(entfac(1),scount(me),
240 & MPI_DOUBLE_PRECISION,entfac(1),
241 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
244 if (me.eq.Master) then
246 write (iout,*) "The FDIMLESS array before sorting"
248 write (iout,*) i,fdimless(i)
255 call mysort1(ntot(islice),Fdimless,iperm)
257 write (iout,*) "The FDIMLESS array after sorting"
259 write (iout,*) i,iperm(i),fdimless(i)
264 qfree=qfree+exp(-fdimless(i)+fdimless(1))
266 c write (iout,*) "qfree",qfree
269 do i=1,min0(ntot(islice),ensembles)
270 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
272 write (iout,*) i,ib,beta_h(ib,iparm),
273 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
274 & potE(iperm(i),iparm),
275 & -entfac(iperm(i)),fdimless(i),sumprob
277 if (sumprob.gt.0.99d0) goto 122
283 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
285 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
290 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
293 if (iproc.ge.nprocs) then
294 write (iout,*) "Fatal error: processor out of range",iproc
299 close (ientout,status="delete")
303 ik=ii-indstart(iproc)+1
304 if (iproc.ne.Master) then
305 if (me.eq.iproc) then
307 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
308 & " energy",potE(ik,iparm)
310 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
311 & Master,i,WHAM_COMM,IERROR)
312 else if (me.eq.Master) then
313 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
314 & WHAM_COMM,STATUS,IERROR)
316 else if (me.eq.Master) then
317 enepot(i)=potE(ik,iparm)
322 enepot(i)=potE(iperm(i),iparm)
326 if (me.eq.Master) then
328 write(licz3,'(bz,i3.3)') iparm
329 write(licz2,'(bz,i2.2)') islice
330 if (temper.lt.100.0d0) then
331 write(ctemper,'(f3.0)') temper
332 else if (temper.lt.1000.0) then
333 write (ctemper,'(f4.0)') temper
335 write (ctemper,'(f5.0)') temper
337 if (nparmset.eq.1) then
338 if (separate_parset) then
339 write(licz4,'(bz,i3.3)') myparm
340 pdbname=prefix(:ilen(prefix))//"_par"//licz4
342 pdbname=prefix(:ilen(prefix))
345 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
347 if (nslice.eq.1) then
348 pdbname=pdbname(:ilen(pdbname))//"_T_"//
349 & ctemper(:ilen(ctemper))//"pdb"
351 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
352 & ctemper(:ilen(ctemper))//"pdb"
354 open(ipdb,file=pdbname)
356 read (ientout,rec=iperm(i))
357 & ((csingle(l,k),l=1,3),k=1,nres),
358 & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
359 & nss,(ihpb(k),jhpb(k),k=1,nss),
360 & eini,efree,rmsdev,iscor
367 call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
377 close(ientout,status="delete")
381 !--------------------------------------------------
382 subroutine mysort1(n, x, ipermut)
384 integer i,j,imax,ipm,n
392 if (x(j).lt.xtemp) then
400 ipermut(imax)=ipermut(i)