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"
14 include "COMMON.HOMOLOGY"
16 include "COMMON.ENERGIES"
17 include "COMMON.FFIELD"
18 include "COMMON.INTERACT"
19 include "COMMON.SBRIDGE"
20 include "COMMON.CHAIN"
21 include "COMMON.PROTFILES"
23 real*4 csingle(3,maxres2)
24 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
25 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
26 double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
27 & escloc,eliptran,esaxs,
28 & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
29 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt,
30 & ehomology_constr,edfadis,edfator,edfanei,edfabet
31 integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
32 double precision qfree,sumprob,eini,efree,rmsdev
34 character*2 licz1,licz2
35 character*3 licz3,licz4
39 real*4 Fdimless(MaxStr),Fdimless_(MaxStr)
40 double precision enepot(MaxStr)
45 if (me.eq.Master) then
47 write (licz2,'(bz,i2.2)') islice
49 if (.not.separate_parset) then
50 bxname = prefix(:ilen(prefix))//".bx"
52 write (licz3,'(bz,i3.3)') myparm
53 bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
56 if (.not.separate_parset) then
57 bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
59 write (licz3,'(bz,i3.3)') myparm
60 bxname = prefix(:ilen(prefix))//"par_"//licz3//
61 & "_slice_"//licz2//".bx"
64 open (ientout,file=bxname,status="unknown",
65 & form="unformatted",access="direct",recl=lenrec1)
70 if (iparm.ne.iparmprint) exit
71 call restore_parm(iparm)
74 write (iout,*) "iparm",iparm," ib",ib
76 temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
77 c quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
84 c fT(l)=kfacl/(kfacl-1.0d0+quotl)
86 if (rescale_mode.eq.1) then
87 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
89 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
90 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
102 fT(l)=kfacl/(kfacl-1.0d0+quotl)
104 else if (rescale_mode.eq.2) then
105 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
107 tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
108 ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
117 fT(l)=1.12692801104297249644d0/
118 & dlog(dexp(quotl)+dexp(-quotl))
120 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
121 else if (rescale_mode.eq.0) then
127 & "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
136 evdw=enetb(1,i,iparm)
137 evdw_t=enetb(21,i,iparm)
139 evdw2_14=enetb(17,i,iparm)
140 evdw2=enetb(2,i,iparm)+evdw2_14
142 evdw2=enetb(2,i,iparm)
147 evdw1=enetb(16,i,iparm)
152 ecorr=enetb(4,i,iparm)
153 ecorr5=enetb(5,i,iparm)
154 ecorr6=enetb(6,i,iparm)
155 eel_loc=enetb(7,i,iparm)
156 eello_turn3=enetb(8,i,iparm)
157 eello_turn4=enetb(9,i,iparm)
158 eturn6=enetb(10,i,iparm)
159 ebe=enetb(11,i,iparm)
160 escloc=enetb(12,i,iparm)
161 etors=enetb(13,i,iparm)
162 etors_d=enetb(14,i,iparm)
163 ehpb=enetb(15,i,iparm)
164 estr=enetb(18,i,iparm)
165 esccor=enetb(19,i,iparm)
166 edihcnstr=enetb(20,i,iparm)
167 eliptran=enetb(22,i,iparm)
168 esaxs=enetb(26,i,iparm)
169 ehomology_constr=enetb(27,i,iparm)
170 edfadis=enetb(28,i,iparm)
171 edfator=enetb(29,i,iparm)
172 edfanei=enetb(30,i,iparm)
173 edfabet=enetb(31,i,iparm)
175 & ehomology_constr=waga_homology(homol_nset)*ehomology_constr
177 if (shield_mode.gt.0) then
178 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
180 & +ft(1)*wvdwpp*evdw1
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+wsaxs*esaxs
190 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
192 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
194 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
195 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
196 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
197 & +ft(2)*wturn3*eello_turn3
198 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
199 & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
200 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
203 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
206 if (shield_mode.gt.0) then
207 etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
208 & +ft(1)*welec*(ees+evdw1)
209 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
210 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
211 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
212 & +ft(2)*wturn3*eello_turn3
213 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
214 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
215 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
218 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
220 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
221 & +ft(1)*welec*(ees+evdw1)
222 & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
223 & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
224 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
225 & +ft(2)*wturn3*eello_turn3
226 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
227 & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
228 & +wbond*estr+wliptran*eliptran+wsaxs*esaxs
231 & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
236 & beta_h(ib,iparm)*etot-entfac(i)
239 write (iout,*) i,indstart(me)+i-1,ib,
240 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
241 & -entfac(i),Fdimless(i)
244 Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
249 call MPI_Gatherv(Fdimless_(1),scount(me),
250 & MPI_REAL,Fdimless(1),
251 & scount(0),idispl(0),MPI_REAL,Master,
254 call MPI_Gatherv(potE(1,iparm),scount(me),
255 & MPI_DOUBLE_PRECISION,potE(1,iparm),
256 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
258 call MPI_Gatherv(entfac(1),scount(me),
259 & MPI_DOUBLE_PRECISION,entfac(1),
260 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
263 if (me.eq.Master) then
265 write (iout,*) "The FDIMLESS array before sorting"
267 write (iout,*) i,fdimless(i)
274 call mysort1(ntot(islice),Fdimless,iperm)
276 write (iout,*) "The FDIMLESS array after sorting"
278 write (iout,*) i,iperm(i),fdimless(i)
283 qfree=qfree+exp(-fdimless(i)+fdimless(1))
285 c write (iout,*) "qfree",qfree
288 do i=1,min0(ntot(islice),ensembles)
289 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
291 write (iout,*) i,ib,beta_h(ib,iparm),
292 & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
293 & potE(iperm(i),iparm),
294 & -entfac(iperm(i)),fdimless(i),sumprob
296 if (sumprob.gt.0.99d0) goto 122
302 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
304 call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
309 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
312 if (iproc.ge.nprocs) then
313 write (iout,*) "Fatal error: processor out of range",iproc
318 close (ientout,status="delete")
322 ik=ii-indstart(iproc)+1
323 if (iproc.ne.Master) then
324 if (me.eq.iproc) then
326 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
327 & " energy",potE(ik,iparm)
329 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
330 & Master,i,WHAM_COMM,IERROR)
331 else if (me.eq.Master) then
332 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
333 & WHAM_COMM,STATUS,IERROR)
335 else if (me.eq.Master) then
336 enepot(i)=potE(ik,iparm)
341 enepot(i)=potE(iperm(i),iparm)
345 if (me.eq.Master) then
347 write(licz3,'(bz,i3.3)') iparm
348 write(licz2,'(bz,i2.2)') islice
349 if (temper.lt.100.0d0) then
350 write(ctemper,'(f3.0)') temper
351 else if (temper.lt.1000.0) then
352 write (ctemper,'(f4.0)') temper
354 write (ctemper,'(f5.0)') temper
356 if (nparmset.eq.1) then
357 if (separate_parset) then
358 write(licz4,'(bz,i3.3)') myparm
359 pdbname=prefix(:ilen(prefix))//"_par"//licz4
361 pdbname=prefix(:ilen(prefix))
364 pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
366 if (nslice.eq.1) then
367 pdbname=pdbname(:ilen(pdbname))//"_T_"//
368 & ctemper(:ilen(ctemper))//"pdb"
370 pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
371 & ctemper(:ilen(ctemper))//"pdb"
373 open(ipdb,file=pdbname)
374 write (iout,*) "Before reading nlist",nlist
376 read (ientout,rec=iperm(i))
377 & ((csingle(l,k),l=1,3),k=1,nres),
378 & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
379 & nss,(ihpb(k),jhpb(k),k=1,nss),
380 & eini,efree,rmsdev,iscor
387 call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
397 close(ientout,status="delete")
401 !--------------------------------------------------
402 subroutine mysort1(n, x, ipermut)
404 integer i,j,imax,ipm,n
412 if (x(j).lt.xtemp) then
420 ipermut(imax)=ipermut(i)