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,
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)
167 if (shield_mode.gt.0) then
168 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
170 & +ft(1)*wvdwpp*evdw1
171 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
172 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
173 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
174 & +ft(2)*wturn3*eello_turn3
175 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
176 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
177 & +wbond*estr+wliptran*eliptran
179 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
181 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
182 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
183 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
184 & +ft(2)*wturn3*eello_turn3
185 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
186 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
187 & +wbond*estr+wliptran*eliptran
190 if (shield_mode.gt.0) then
191 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
192 & +ft(1)*welec*(ees+evdw1)
193 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
194 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
195 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
196 & +ft(2)*wturn3*eello_turn3
197 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
198 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
199 & +wbond*estr+wliptran*eliptran
201 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
202 & +ft(1)*welec*(ees+evdw1)
203 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
204 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
205 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
206 & +ft(2)*wturn3*eello_turn3
207 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
208 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
209 & +wbond*estr+wliptran*eliptran
215 & beta_h(ib,iparm)*etot-entfac(i)
218 write (iout,*) i,indstart(me)+i-1,ib,
219 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
220 & -entfac(i),Fdimless(i)
223 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
228 call MPI_Gatherv(Fdimless(1),scount(me),
229 & MPI_REAL,Fdimless(1),
230 & scount(0),idispl(0),MPI_REAL,Master,
233 call MPI_Gatherv(potE(1,iparm),scount(me),
234 & MPI_DOUBLE_PRECISION,potE(1,iparm),
235 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
237 call MPI_Gatherv(entfac(1),scount(me),
238 & MPI_DOUBLE_PRECISION,entfac(1),
239 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
242 if (me.eq.Master) then
244 write (iout,*) "The FDIMLESS array before sorting"
246 write (iout,*) i,fdimless(i)
253 call mysort1(ntot(islice),Fdimless,iperm)
255 write (iout,*) "The FDIMLESS array after sorting"
257 write (iout,*) i,iperm(i),fdimless(i)
262 qfree=qfree+exp(-fdimless(i)+fdimless(1))
264 c write (iout,*) "qfree",qfree
267 do i=1,min0(ntot(islice),ensembles)
268 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
270 write (iout,*) i,ib,beta_h(ib,iparm),
271 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
272 & potE(iperm(i),iparm),
273 & -entfac(iperm(i)),fdimless(i),sumprob
275 if (sumprob.gt.0.99d0) goto 122
281 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
283 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
288 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
291 if (iproc.ge.nprocs) then
292 write (iout,*) "Fatal error: processor out of range",iproc
297 close (ientout,status="delete")
301 ik=ii-indstart(iproc)+1
302 if (iproc.ne.Master) then
303 if (me.eq.iproc) then
305 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
306 & " energy",potE(ik,iparm)
308 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
309 & Master,i,WHAM_COMM,IERROR)
310 else if (me.eq.Master) then
311 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
312 & WHAM_COMM,STATUS,IERROR)
314 else if (me.eq.Master) then
315 enepot(i)=potE(ik,iparm)
320 enepot(i)=potE(iperm(i),iparm)
324 if (me.eq.Master) then
326 write(licz3,'(bz,i3.3)') iparm
327 write(licz2,'(bz,i2.2)') islice
328 if (temper.lt.100.0d0) then
329 write(ctemper,'(f3.0)') temper
330 else if (temper.lt.1000.0) then
331 write (ctemper,'(f4.0)') temper
333 write (ctemper,'(f5.0)') temper
335 if (nparmset.eq.1) then
336 if (separate_parset) then
337 write(licz4,'(bz,i3.3)') myparm
338 pdbname=prefix(:ilen(prefix))//"_par"//licz4
340 pdbname=prefix(:ilen(prefix))
343 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
345 if (nslice.eq.1) then
346 pdbname=pdbname(:ilen(pdbname))//"_T_"//
347 & ctemper(:ilen(ctemper))//"pdb"
349 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
350 & ctemper(:ilen(ctemper))//"pdb"
352 open(ipdb,file=pdbname)
354 read (ientout,rec=iperm(i))
355 & ((csingle(l,k),l=1,3),k=1,nres),
356 & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
357 & nss,(ihpb(k),jhpb(k),k=1,nss),
358 & eini,efree,rmsdev,iscor
365 call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
375 close(ientout,status="delete")
379 !--------------------------------------------------
380 subroutine mysort1(n, x, ipermut)
382 integer i,j,imax,ipm,n
390 if (x(j).lt.xtemp) then
398 ipermut(imax)=ipermut(i)