1 subroutine probabl(ib,nlist,ncon,*)
2 ! construct the conformational ensembles at REMD temperatures
9 integer ierror,errcode,status(MPI_STATUS_SIZE)
11 include "COMMON.IOUNITS"
13 include "COMMON.FFIELD"
14 include "COMMON.INTERACT"
15 include "COMMON.SBRIDGE"
16 include "COMMON.CHAIN"
17 include "COMMON.CLUSTER"
18 real*4 csingle(3,maxres2)
19 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
20 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
21 double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
22 & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
23 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
25 integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
26 double precision qfree,sumprob,eini,efree,rmsdev
32 real*4 Fdimless(maxconf),Fdimless_(maxconf)
33 double precision energia(0:max_ene)
34 double precision totfree_(maxconf),entfac_(maxconf)
39 c write (iout,*) i,list_conf(i)
42 write (iout,*) me," indstart",indstart(me)," indend",indend(me)
43 call daread_ccoords(indstart(me),indend(me))
45 c write (iout,*) "ncon",ncon
46 temper=1.0d0/(beta_h(ib)*1.987D-3)
47 c write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
48 c quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
55 c fT(l)=kfacl/(kfacl-1.0d0+quotl)
57 if (rescale_mode.eq.1) then
58 quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
65 fT(l)=kfacl/(kfacl-1.0d0+quotl)
68 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
70 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
71 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
72 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
83 else if (rescale_mode.eq.2) then
84 quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
88 fT(l)=1.12692801104297249644d0/
89 & dlog(dexp(quotl)+dexp(-quotl))
91 c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
94 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
96 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
97 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
98 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
117 c write (iout,*) "i",i," ii",ii
122 c(k,j)=allcart(k,j,i)
123 c(k,j+nres)=allcart(k,j+nres,i)
128 c(k,nres+nres)=c(k,nres)
132 ihpb(j)=ihpb_all(j,i)
133 jhpb(j)=jhpb_all(j,i)
135 call int_from_cart1(.false.)
136 call etotal(energia(0),fT)
137 totfree(i)=energia(0)
138 c write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
139 c write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
140 c call enerprint(energia(0),fT)
141 c call pdbout(totfree(i),16,i)
143 write (iout,*) i," energia",(energia(j),j=0,max_ene)
144 write (iout,*) "etot", etot
145 write (iout,*) "ft(6)", ft(6)
148 enetb(k,i)=energia(k)
152 c write (iout,*) evdw
156 evdw2=enetb(2,i)+evdw2_14
172 eello_turn3=enetb(8,i)
173 eello_turn4=enetb(9,i)
182 edihcnstr=enetb(20,i)
185 c etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
186 c &ft(1)*welec*ees+wvdwpp*evdw1
187 c & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
188 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
189 c & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
190 c & +ft(2)*wturn3*eello_turn3
191 c & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
192 c & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
195 c etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
196 c & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
197 c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
198 c & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
199 c & +ft(2)*wturn3*eello_turn3
200 c & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
201 c & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
205 write (iout,*) "etot2", etot
206 write (iout,*) "evdw", wsc, evdw,evdw_t
207 write (iout,*) "evdw2", wscp, evdw2
208 write (iout,*) "welec", ft(1),welec,ees
209 write (iout,*) "evdw1", wvdwpp,evdw1
210 write (iout,*) "ebe" ebe,wang
213 Fdimless_(i)=beta_h(ib)*etot+entfac(ii)
216 Fdimless(i)=beta_h(ib)*etot+entfac(ii)
220 write (iout,*) "fdim calc", i,ii,ib,
221 & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
222 & entfac(ii),Fdimless(i)
226 call MPI_Gatherv(Fdimless_(1),scount(me),
227 & MPI_REAL,Fdimless(1),
228 & scount(0),idispl(0),MPI_REAL,Master,
229 & MPI_COMM_WORLD, IERROR)
230 call MPI_Gatherv(totfree_(1),scount(me),
231 & MPI_DOUBLE_PRECISION,totfree(1),
232 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
233 & MPI_COMM_WORLD, IERROR)
234 call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
235 & MPI_DOUBLE_PRECISION,entfac_(1),
236 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
237 & MPI_COMM_WORLD, IERROR)
238 if (me.eq.Master) then
244 write (iout,*) "The FDIMLESS array before sorting"
246 c write (iout,*) i,fdimless(i)
249 call mysort1(ncon,Fdimless,list_conf)
251 write (iout,*) "The FDIMLESS array after sorting"
253 write (iout,*) i,list_conf(i),fdimless(i)
257 totfree(i)=fdimless(i)
261 qfree=qfree+exp(-fdimless(i)+fdimless(1))
262 c write (iout,*) "fdimless", fdimless(i)
264 c write (iout,*) "qfree",qfree
267 write (iout,*) "ncon", ncon,maxstr_proc
268 do i=1,min0(ncon,maxstr_proc)-1
269 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
271 write (iout,*) i,ib,beta_h(ib),
272 & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
273 & totfree(list_conf(i)),
274 & -entfac(list_conf(i)),fdimless(i),sumprob
276 if (sumprob.gt.prob_limit) goto 122
277 c if (sumprob.gt.1.00d0) goto 122
283 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
285 call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
288 c write (iout,*) "iproc",iproc," indstart",indstart(iproc),
289 c & " indend",indend(iproc)
291 write (iout,*) "nlist",nlist
295 !--------------------------------------------------
296 subroutine mysort1(n, x, ipermut)
298 integer i,j,imax,ipm,n
306 if (x(j).lt.xtemp) then
314 ipermut(imax)=ipermut(i)