543422a35f71d8c3df6be90f043ad1134e253206
[unres.git] / probabl.F
1       subroutine probabl(ib,nlist,ncon,*)
2 ! construct the conformational ensembles at REMD temperatures
3       implicit none
4       include "DIMENSIONS"
5       include "sizesclu.dat"
6 #ifdef MPI
7       include "mpif.h"
8       include "COMMON.MPI"
9       integer ierror,errcode,status(MPI_STATUS_SIZE) 
10 #endif
11       include "COMMON.IOUNITS"
12       include "COMMON.FREE"
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,
24      &      evdw_t
25       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
26       double precision qfree,sumprob,eini,efree,rmsdev
27       character*80 bxname
28       character*2 licz1
29       character*5 ctemper
30       integer ilen,ijk
31       external ilen
32       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
33       double precision energia(0:max_ene), totfree_buf(0:maxconf),
34      &  entfac_buf(maxconf)
35       do i=1,ncon
36         list_conf(i)=i
37       enddo
38 c      do i=1,ncon
39 c        write (iout,*) i,list_conf(i)
40 c      enddo
41 #ifdef MPI
42       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
43       call daread_ccoords(indstart(me),indend(me))
44 #endif
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)
49 c      quotl=1.0d0
50 c      kfacl=1.0d0
51 c      do l=1,5
52 c        quotl1=quotl
53 c        quotl=quotl*quot
54 c        kfacl=kfacl*kfac
55 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
56 c      enddo
57             if (rescale_mode.eq.1) then
58               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
59               quotl=1.0d0
60               kfacl=1.0d0
61               do l=1,5
62                 quotl1=quotl
63                 quotl=quotl*quot
64                 kfacl=kfacl*kfac
65                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
66               enddo
67 #if defined(FUNCTH)
68               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
69      &                  320.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)
73 #elif defined(FUNCT)
74               fT(6)=betaT/T0
75               ftprim(6)=1.0d0/T0
76               ftbis(6)=0.0d0
77 #else
78               fT(6)=1.0d0
79               ftprim(6)=0.0d0
80               ftbis(6)=0.0d0
81 #endif
82
83             else if (rescale_mode.eq.2) then
84               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
85               quotl=1.0d0
86               do l=1,5
87                 quotl=quotl*quot
88                 fT(l)=1.12692801104297249644d0/
89      &             dlog(dexp(quotl)+dexp(-quotl))
90               enddo
91 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
92 c              call flush(iout)
93 #if defined(FUNCTH)
94               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
95      &                  320.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)
99 #elif defined(FUNCT)
100               fT(6)=betaT/T0
101               ftprim(6)=1.0d0/T0
102               ftbis(6)=0.0d0
103 #else
104               fT(6)=1.0d0
105               ftprim(6)=0.0d0
106               ftbis(6)=0.0d0
107 #endif
108             endif
109
110 #ifdef MPI
111       do i=1,scount(me)
112         ii=i+indstart(me)-1
113 #else
114       do i=1,ncon
115         ii=i
116 #endif
117 c        write (iout,*) "i",i," ii",ii
118 c        call flush(iout)
119         if (ib.eq.1) then
120           do j=1,nres
121             do k=1,3
122               c(k,j)=allcart(k,j,i)
123               c(k,j+nres)=allcart(k,j+nres,i)
124             enddo
125           enddo
126           do k=1,3
127             c(k,nres+1)=c(k,1)
128             c(k,nres+nres)=c(k,nres)
129           enddo
130           nss=nss_all(i)
131           do j=1,nss
132             ihpb(j)=ihpb_all(j,i)
133             jhpb(j)=jhpb_all(j,i)
134           enddo 
135           call int_from_cart1(.false.)
136           call etotal(energia(0),fT)
137           totfree(i)=energia(0)         
138           totfree_buf(i)=totfree(i)
139 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
140 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
141 c          call enerprint(energia(0),fT)
142 c          call pdbout(totfree(i),16,i)
143 #ifdef DEBUG
144           write (iout,*) i," energia",(energia(j),j=0,19)
145           write (iout,*) "etot", etot
146           write (iout,*) "ft(6)", ft(6)
147 #endif
148           do k=1,max_ene
149             enetb(k,i)=energia(k)
150           enddo
151         endif
152         evdw=enetb(1,i)
153 c        write (iout,*) evdw
154        etot=energia(0)
155 #ifdef SCP14
156         evdw2_14=enetb(17,i)
157         evdw2=enetb(2,i)+evdw2_14
158 #else
159         evdw2=enetb(2,i)
160         evdw2_14=0.0d0
161 #endif
162 #ifdef SPLITELE
163         ees=enetb(3,i)
164         evdw1=enetb(16,i)
165 #else
166         ees=enetb(3,i)
167         evdw1=0.0d0
168 #endif
169         ecorr=enetb(4,i)
170         ecorr5=enetb(5,i)
171         ecorr6=enetb(6,i)
172         eel_loc=enetb(7,i)
173         eello_turn3=enetb(8,i)
174         eello_turn4=enetb(9,i)
175         eturn6=enetb(10,i)
176         ebe=enetb(11,i)
177         escloc=enetb(12,i)
178         etors=enetb(13,i)
179         etors_d=enetb(14,i)
180         ehpb=enetb(15,i)
181         estr=enetb(18,i)
182         esccor=enetb(19,i)
183         edihcnstr=enetb(20,i)
184 c#ifdef SPLITELE
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
193 c     &  +wbond*estr
194 c#else
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
202 c     &  +wbond*estr
203 c#endif
204 #ifdef DEBUG
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
211 #endif        
212         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
213         Fdimless_buf(i)=Fdimless(i)
214         totfree(i)=etot
215         totfree_buf(i)=totfree(i)
216 #ifdef DEBUG
217         write (iout,*) "fdim calc", i,ii,ib,
218      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
219      &   entfac(ii),Fdimless(i)
220 #endif
221       enddo   ! i
222
223       do ijk=1,maxconf
224       entfac_buf(ijk)=entfac(ijk)
225       Fdimless_buf(ijk)=Fdimless(ijk)
226       enddo
227       do ijk=0,maxconf
228       totfree_buf(ijk)=totfree(ijk)
229       enddo
230
231
232 c      scount_buf=scount(me)
233 c      scount_buf2=scount(0)
234
235 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
236
237 #ifdef MPI
238       WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
239       call MPI_Gatherv(Fdimless_buf(1),scount(me),
240      & MPI_REAL,Fdimless(1),
241      & scount(0),idispl(0),MPI_REAL,Master,
242      & MPI_COMM_WORLD, IERROR)
243       WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
244       call MPI_Gatherv(totfree_buf(1),scount(me),
245      & MPI_DOUBLE_PRECISION,totfree(1),
246      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
247      & MPI_COMM_WORLD, IERROR)
248       WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
249       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
250      & MPI_DOUBLE_PRECISION,entfac(1),
251      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
252      & MPI_COMM_WORLD, IERROR)
253       WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
254       if (me.eq.Master) then
255       WRITE (iout,*) "me.eq.Master"
256 #endif
257 #ifdef DEBUG
258         write (iout,*) "The FDIMLESS array before sorting"
259         do i=1,ncon
260 c          write (iout,*) i,fdimless(i)
261         enddo
262 #endif
263       WRITE (iout,*) "Wchodze do call mysort1"
264         call mysort1(ncon,Fdimless,list_conf)
265       WRITE (iout,*) "Wychodze z call mysort1"
266 #ifdef DEBUG
267         write (iout,*) "The FDIMLESS array after sorting"
268         do i=1,ncon
269           write (iout,*) i,list_conf(i),fdimless(i)
270         enddo
271 #endif
272       WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
273         do i=1,ncon
274           totfree(i)=fdimless(i)
275         enddo
276         qfree=0.0d0
277         do i=1,ncon
278           qfree=qfree+exp(-fdimless(i)+fdimless(1))
279 c          write (iout,*) "fdimless", fdimless(i)
280         enddo
281 c        write (iout,*) "qfree",qfree
282         nlist=1
283         sumprob=0.0
284         write (iout,*) "ncon", ncon,maxstr_proc
285         do i=1,min0(ncon,maxstr_proc)-1 
286           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
287 #ifdef DEBUG
288           write (iout,*) i,ib,beta_h(ib),
289      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
290      &     totfree(list_conf(i)),
291      &     -entfac(list_conf(i)),fdimless(i),sumprob
292 #endif
293           if (sumprob.gt.prob_limit) goto 122
294 c          if (sumprob.gt.1.00d0) goto 122
295           nlist=nlist+1
296         enddo  
297   122   continue
298 #ifdef MPI
299       endif
300       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, 
301      &   IERROR)
302       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
303      &   IERROR)
304 c      do iproc=0,nprocs
305 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
306 c     &   " indend",indend(iproc) 
307 c      enddo
308       write (iout,*) "nlist",nlist
309 #endif
310       return
311       end
312 !--------------------------------------------------
313       subroutine mysort1(n, x, ipermut)
314       implicit none
315       integer i,j,imax,ipm,n
316       real x(n)
317       integer ipermut(n)
318       real xtemp
319       do i=1,n
320         xtemp=x(i)
321         imax=i
322         do j=i+1,n
323           if (x(j).lt.xtemp) then
324             imax=j
325             xtemp=x(j)
326           endif
327         enddo
328         x(imax)=x(i)
329         x(i)=xtemp
330         ipm=ipermut(imax)
331         ipermut(imax)=ipermut(i)
332         ipermut(i)=ipm
333       enddo
334       return
335       end