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.CONTROL"
12 include "COMMON.IOUNITS"
14 include "COMMON.FFIELD"
15 include "COMMON.INTERACT"
16 include "COMMON.SBRIDGE"
17 include "COMMON.CHAIN"
18 include "COMMON.CLUSTER"
19 real*4 csingle(3,maxres2)
20 double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
21 & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
22 double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
23 & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
24 & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
25 & evdw_t,esaxs,eliptran,ethetacnstr,ehomology_constr,
26 & edfadis,edfator,edfanei,edfabet
27 integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
28 double precision qfree,sumprob,eini,efree,rmsdev
34 character*80 structure/'Structure'/
35 real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
36 double precision energia(0:max_ene), totfree_buf(0:maxconf),
38 double precision buffer(maxconf)
43 c write (iout,*) i,list_conf(i)
46 write (iout,*) me," indstart",indstart(me)," indend",indend(me)
47 call daread_ccoords(indstart(me),indend(me))
49 C write (iout,*) "ncon",ncon
51 temper=1.0d0/(beta_h(ib)*1.987D-3)
52 if (rescale_mode.eq.1) then
53 quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
60 fT(l)=kfacl/(kfacl-1.0d0+quotl)
63 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
65 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
66 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
67 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
78 else if (rescale_mode.eq.2) then
79 quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
83 fT(l)=1.12692801104297249644d0/
84 & dlog(dexp(quotl)+dexp(-quotl))
86 c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
89 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
91 ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
92 ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
93 & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
112 C write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
117 c(k,j)=allcart(k,j,i)
118 c(k,j+nres)=allcart(k,j+nres,i)
119 C write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
120 C & c(k,j+nres),allcart(k,j+nres,i)
123 C write(iout,*) "out of j loop"
127 c(k,nres+nres)=c(k,nres)
129 C write(iout,*) "after nres+nres",nss_all(i)
133 ihpb(j)=ihpb_all(j,i)
134 jhpb(j)=jhpb_all(j,i)
136 call int_from_cart1(.false.)
137 call etotal(energia(0),fT)
139 write (structure(9:),'(bz,i6.6)') i
140 call TMscore_sub(rmsdev,gdt_ts_tb(i),
141 & gdt_ha_tb(i),tmscore_tb(i),Structure,.false.)
143 write (iout,*) rmsdev,gdt_ts_tb(i),gdt_ha_tb(i),
147 totfree(i)=energia(0)
148 totfree_buf(i)=totfree(i)
149 c write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
150 c write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
151 c call pdbout(totfree(i),16,i)
154 write (iout,*) "conformation", i
155 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
156 & ((c(l,k+nres),l=1,3),k=nnt,nct)
157 call enerprint(energia(0),fT)
160 Fdimless(i)=beta_h(ib)*etot+entfac(ii)
161 Fdimless_buf(i)=Fdimless(i)
163 totfree_buf(i)=totfree(i)
165 write (iout,*) "fdim calc", i,ii,ib,
166 & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
167 & entfac(ii),Fdimless(i)
172 entfac_buf(ijk)=entfac(ijk)
173 Fdimless_buf(ijk)=Fdimless(ijk)
176 totfree_buf(ijk)=totfree(ijk)
180 c scount_buf=scount(me)
181 c scount_buf2=scount(0)
183 c entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
186 c WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
187 call MPI_Gatherv(Fdimless_buf(1),scount(me),
188 & MPI_REAL,Fdimless(1),
189 & scount(0),idispl(0),MPI_REAL,Master,
190 & MPI_COMM_WORLD, IERROR)
191 c WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
192 call MPI_Gatherv(totfree_buf(1),scount(me),
193 & MPI_DOUBLE_PRECISION,totfree(1),
194 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
195 & MPI_COMM_WORLD, IERROR)
196 c WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
197 call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
198 & MPI_DOUBLE_PRECISION,entfac(1),
199 & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
200 & MPI_COMM_WORLD, IERROR)
201 c WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
204 buffer(i)=gdt_ts_tb(i)
206 call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
207 & gdt_ts_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
208 & MPI_COMM_WORLD,IERROR)
210 buffer(i)=gdt_ha_tb(i)
212 call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
213 & gdt_ha_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
214 & MPI_COMM_WORLD,IERROR)
216 buffer(i)=tmscore_tb(i)
218 call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
219 & tmscore_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
220 & MPI_COMM_WORLD,IERROR)
222 if (me.eq.Master) then
223 c WRITE (iout,*) "me.eq.Master"
226 write (iout,*) "The FDIMLESS array before sorting"
228 write (iout,*) i,fdimless(i)
231 c WRITE (iout,*) "Wchodze do call mysort1"
232 call mysort1(ncon,Fdimless,list_conf)
233 c WRITE (iout,*) "Wychodze z call mysort1"
235 write (iout,*) "The FDIMLESS array after sorting"
237 write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
238 & gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
241 c WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
243 totfree(i)=fdimless(i)
247 qfree=qfree+exp(-fdimless(i)+fdimless(1))
248 c write (iout,*) "fdimless", fdimless(i)
250 c write (iout,*) "qfree",qfree
253 write (iout,*) "ncon", ncon,maxstr_proc
254 do i=1,min0(ncon,maxstr_proc)-1
255 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
257 write (iout,*) i,ib,beta_h(ib),
258 & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
259 & totfree(list_conf(i)),
260 & -entfac(list_conf(i)),fdimless(i),sumprob
262 if (sumprob.gt.prob_limit) goto 122
263 c if (sumprob.gt.1.00d0) goto 122
269 call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
271 call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
274 c write (iout,*) "iproc",iproc," indstart",indstart(iproc),
275 c & " indend",indend(iproc)
277 write (iout,*) "nlist",nlist
281 !--------------------------------------------------
282 subroutine mysort1(n, x, ipermut)
284 integer i,j,imax,ipm,n
292 if (x(j).lt.xtemp) then
300 ipermut(imax)=ipermut(i)