rm uscdiff
[unres.git] / source / cluster / wham / src-M / 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 C      call flush(iout)
47       temper=1.0d0/(beta_h(ib)*1.987D-3)
48 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
49 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
50 c      quotl=1.0d0
51 c      kfacl=1.0d0
52 c      do l=1,5
53 c        quotl1=quotl
54 c        quotl=quotl*quot
55 c        kfacl=kfacl*kfac
56 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
57 c      enddo
58 C#define DEBUG
59             if (rescale_mode.eq.1) then
60               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
61               quotl=1.0d0
62               kfacl=1.0d0
63               do l=1,5
64                 quotl1=quotl
65                 quotl=quotl*quot
66                 kfacl=kfacl*kfac
67                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
68               enddo
69 #if defined(FUNCTH)
70               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
71      &                  320.0d0
72               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
73              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
74      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
75 #elif defined(FUNCT)
76               fT(6)=betaT/T0
77               ftprim(6)=1.0d0/T0
78               ftbis(6)=0.0d0
79 #else
80               fT(6)=1.0d0
81               ftprim(6)=0.0d0
82               ftbis(6)=0.0d0
83 #endif
84
85             else if (rescale_mode.eq.2) then
86               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
87               quotl=1.0d0
88               do l=1,5
89                 quotl=quotl*quot
90                 fT(l)=1.12692801104297249644d0/
91      &             dlog(dexp(quotl)+dexp(-quotl))
92               enddo
93 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
94 c              call flush(iout)
95 #if defined(FUNCTH)
96               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
97      &                  320.0d0
98               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
99              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
100      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
101 #elif defined(FUNCT)
102               fT(6)=betaT/T0
103               ftprim(6)=1.0d0/T0
104               ftbis(6)=0.0d0
105 #else
106               fT(6)=1.0d0
107               ftprim(6)=0.0d0
108               ftbis(6)=0.0d0
109 #endif
110             endif
111
112 #ifdef MPI
113       do i=1,scount(me)
114         ii=i+indstart(me)-1
115 #else
116       do i=1,ncon
117         ii=i
118 #endif
119 C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
120 c        call flush(iout)
121         if (ib.eq.1) then
122           do j=1,nres
123             do k=1,3
124               c(k,j)=allcart(k,j,i)
125               c(k,j+nres)=allcart(k,j+nres,i)
126 C              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
127 C     &        c(k,j+nres),allcart(k,j+nres,i)
128             enddo
129           enddo
130 C          write(iout,*) "out of j loop"
131 C          call flush(iout)
132           do k=1,3
133             c(k,nres+1)=c(k,1)
134             c(k,nres+nres)=c(k,nres)
135           enddo
136 C          write(iout,*) "after nres+nres",nss_all(i)
137 C          call flush(iout)
138           nss=nss_all(i)
139           do j=1,nss
140             ihpb(j)=ihpb_all(j,i)
141             jhpb(j)=jhpb_all(j,i)
142           enddo 
143           call int_from_cart1(.false.)
144 C          write(iout,*) "before etotal"
145 C          call flush(iout)
146           call etotal(energia(0),fT)
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 enerprint(energia(0),fT)
152 c          call pdbout(totfree(i),16,i)
153 #ifdef DEBUG
154           write (iout,*) i," energia",(energia(j),j=0,max_ene)
155           write (iout,*) "etot", etot
156           write (iout,*) "ft(6)", ft(6)
157 #endif
158           do k=1,max_ene
159             enetb(k,i)=energia(k)
160           enddo
161         endif
162         evdw=enetb(1,i)
163 c        write (iout,*) evdw
164        etot=energia(0)
165 #ifdef SCP14
166         evdw2_14=enetb(17,i)
167         evdw2=enetb(2,i)+evdw2_14
168 #else
169         evdw2=enetb(2,i)
170         evdw2_14=0.0d0
171 #endif
172 #ifdef SPLITELE
173         ees=enetb(3,i)
174         evdw1=enetb(16,i)
175 #else
176         ees=enetb(3,i)
177         evdw1=0.0d0
178 #endif
179         ecorr=enetb(4,i)
180         ecorr5=enetb(5,i)
181         ecorr6=enetb(6,i)
182         eel_loc=enetb(7,i)
183         eello_turn3=enetb(8,i)
184         eello_turn4=enetb(9,i)
185         eturn6=enetb(10,i)
186         ebe=enetb(11,i)
187         escloc=enetb(12,i)
188         etors=enetb(13,i)
189         etors_d=enetb(14,i)
190         ehpb=enetb(15,i)
191         estr=enetb(18,i)
192         esccor=enetb(19,i)
193         edihcnstr=enetb(20,i)
194 c#ifdef SPLITELE
195 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
196 c     &ft(1)*welec*ees+wvdwpp*evdw1
197 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
198 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
199 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
200 c     &  +ft(2)*wturn3*eello_turn3
201 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
202 c     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
203 c     &  +wbond*estr
204 c#else
205 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
206 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
207 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
208 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
209 c     &  +ft(2)*wturn3*eello_turn3
210 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
211 c     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
212 c     &  +wbond*estr
213 c#endif
214 #ifdef DEBUG
215         write (iout,*) "etot2", etot
216         write (iout,*) "evdw", wsc, evdw,evdw_t
217         write (iout,*) "evdw2", wscp, evdw2
218         write (iout,*) "welec", ft(1),welec,ees
219         write (iout,*) "evdw1", wvdwpp,evdw1
220         write (iout,*) "ebe", ebe,wang
221 #endif        
222         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
223         Fdimless_buf(i)=Fdimless(i)
224         totfree(i)=etot
225         totfree_buf(i)=totfree(i)
226 #ifdef DEBUG
227         write (iout,*) "fdim calc", i,ii,ib,
228      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
229      &   entfac(ii),Fdimless(i)
230 #endif
231       enddo   ! i
232
233       do ijk=1,maxconf
234       entfac_buf(ijk)=entfac(ijk)
235       Fdimless_buf(ijk)=Fdimless(ijk)
236       enddo
237       do ijk=0,maxconf
238       totfree_buf(ijk)=totfree(ijk)
239       enddo
240
241
242 c      scount_buf=scount(me)
243 c      scount_buf2=scount(0)
244
245 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
246
247 #ifdef MPI
248 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
249       call MPI_Gatherv(Fdimless_buf(1),scount(me),
250      & MPI_REAL,Fdimless(1),
251      & scount(0),idispl(0),MPI_REAL,Master,
252      & MPI_COMM_WORLD, IERROR)
253 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
254       call MPI_Gatherv(totfree_buf(1),scount(me),
255      & MPI_DOUBLE_PRECISION,totfree(1),
256      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
257      & MPI_COMM_WORLD, IERROR)
258 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
259       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
260      & MPI_DOUBLE_PRECISION,entfac(1),
261      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
262      & MPI_COMM_WORLD, IERROR)
263 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
264       if (me.eq.Master) then
265 c      WRITE (iout,*) "me.eq.Master"
266 #endif
267 #ifdef DEBUG
268         write (iout,*) "The FDIMLESS array before sorting"
269         do i=1,ncon
270 c          write (iout,*) i,fdimless(i)
271         enddo
272 #endif
273 c      WRITE (iout,*) "Wchodze do call mysort1"
274         call mysort1(ncon,Fdimless,list_conf)
275 c      WRITE (iout,*) "Wychodze z call mysort1"
276 #ifdef DEBUG
277         write (iout,*) "The FDIMLESS array after sorting"
278         do i=1,ncon
279           write (iout,*) i,list_conf(i),fdimless(i)
280         enddo
281 #endif
282 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
283         do i=1,ncon
284           totfree(i)=fdimless(i)
285         enddo
286         qfree=0.0d0
287         do i=1,ncon
288           qfree=qfree+exp(-fdimless(i)+fdimless(1))
289 c          write (iout,*) "fdimless", fdimless(i)
290         enddo
291 c        write (iout,*) "qfree",qfree
292         nlist=1
293         sumprob=0.0
294         write (iout,*) "ncon", ncon,maxstr_proc
295         do i=1,min0(ncon,maxstr_proc)-1 
296           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
297 #ifdef DEBUG
298           write (iout,*) "tu szukaj ponizej 7"
299           write (iout,*) i,ib,beta_h(ib),
300      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
301      &     totfree(list_conf(i)),
302      &     -entfac(list_conf(i)),fdimless(i),sumprob
303 #endif
304           if (sumprob.gt.prob_limit) goto 122
305 c          if (sumprob.gt.1.00d0) goto 122
306           nlist=nlist+1
307         enddo  
308   122   continue
309 #ifdef MPI
310       endif
311       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, 
312      &   IERROR)
313       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
314      &   IERROR)
315 c      do iproc=0,nprocs
316 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
317 c     &   " indend",indend(iproc) 
318 c      enddo
319       write (iout,*) "nlist",nlist
320 #endif
321       return
322       end
323 !--------------------------------------------------
324       subroutine mysort1(n, x, ipermut)
325       implicit none
326       integer i,j,imax,ipm,n
327       real x(n)
328       integer ipermut(n)
329       real xtemp
330       do i=1,n
331         xtemp=x(i)
332         imax=i
333         do j=i+1,n
334           if (x(j).lt.xtemp) then
335             imax=j
336             xtemp=x(j)
337           endif
338         enddo
339         x(imax)=x(i)
340         x(i)=xtemp
341         ipm=ipermut(imax)
342         ipermut(imax)=ipermut(i)
343         ipermut(i)=ipm
344       enddo
345       return
346       end