correction to the last commit - cleaning of homo restraints
[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               write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
94               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         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           call enerprint(energia(0),fT)
152 c          call pdbout(totfree(i),16,i)
153 #define DEBUG
154 #ifdef DEBUG
155           write (iout,*) i," energia",(energia(j),j=0,max_ene)
156           write (iout,*) "etot", etot
157           write (iout,*) "ft(6)", ft(6)
158 #endif
159 #undef DEBUG
160           do k=1,max_ene
161             enetb(k,i)=energia(k)
162           enddo
163         endif
164         evdw=enetb(1,i)
165 c        write (iout,*) evdw
166        etot=energia(0)
167 #ifdef SCP14
168         evdw2_14=enetb(17,i)
169         evdw2=enetb(2,i)+evdw2_14
170 #else
171         evdw2=enetb(2,i)
172         evdw2_14=0.0d0
173 #endif
174 #ifdef SPLITELE
175         ees=enetb(3,i)
176         evdw1=enetb(16,i)
177 #else
178         ees=enetb(3,i)
179         evdw1=0.0d0
180 #endif
181         ecorr=enetb(4,i)
182         ecorr5=enetb(5,i)
183         ecorr6=enetb(6,i)
184         eel_loc=enetb(7,i)
185         eello_turn3=enetb(8,i)
186         eello_turn4=enetb(9,i)
187         eturn6=enetb(10,i)
188         ebe=enetb(11,i)
189         escloc=enetb(12,i)
190         etors=enetb(13,i)
191         etors_d=enetb(14,i)
192         ehpb=enetb(15,i)
193         estr=enetb(18,i)
194         esccor=enetb(19,i)
195         edihcnstr=enetb(20,i)
196 c#ifdef SPLITELE
197 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
198 c     &ft(1)*welec*ees+wvdwpp*evdw1
199 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
200 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
201 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
202 c     &  +ft(2)*wturn3*eello_turn3
203 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
204 c     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
205 c     &  +wbond*estr
206 c#else
207 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
208 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
209 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
210 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
211 c     &  +ft(2)*wturn3*eello_turn3
212 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
213 c     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
214 c     &  +wbond*estr
215 c#endif
216 #ifdef DEBUG
217         write (iout,*) "etot2", etot
218         write (iout,*) "evdw", wsc, evdw,evdw_t
219         write (iout,*) "evdw2", wscp, evdw2
220         write (iout,*) "welec", ft(1),welec,ees
221         write (iout,*) "evdw1", wvdwpp,evdw1
222         write (iout,*) "ebe", ebe,wang
223 #endif        
224         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
225         Fdimless_buf(i)=Fdimless(i)
226         totfree(i)=etot
227         totfree_buf(i)=totfree(i)
228 #ifdef DEBUG
229         write (iout,*) "fdim calc", i,ii,ib,
230      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
231      &   entfac(ii),Fdimless(i)
232 #endif
233       enddo   ! i
234
235       do ijk=1,maxconf
236       entfac_buf(ijk)=entfac(ijk)
237       Fdimless_buf(ijk)=Fdimless(ijk)
238       enddo
239       do ijk=0,maxconf
240       totfree_buf(ijk)=totfree(ijk)
241       enddo
242
243
244 c      scount_buf=scount(me)
245 c      scount_buf2=scount(0)
246
247 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
248
249 #ifdef MPI
250       WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
251       call MPI_Gatherv(Fdimless_buf(1),scount(me),
252      & MPI_REAL,Fdimless(1),
253      & scount(0),idispl(0),MPI_REAL,Master,
254      & MPI_COMM_WORLD, IERROR)
255       WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
256       call MPI_Gatherv(totfree_buf(1),scount(me),
257      & MPI_DOUBLE_PRECISION,totfree(1),
258      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
259      & MPI_COMM_WORLD, IERROR)
260       WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
261       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
262      & MPI_DOUBLE_PRECISION,entfac(1),
263      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
264      & MPI_COMM_WORLD, IERROR)
265       WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
266       if (me.eq.Master) then
267       WRITE (iout,*) "me.eq.Master"
268 #endif
269 #ifdef DEBUG
270         write (iout,*) "The FDIMLESS array before sorting"
271         do i=1,ncon
272 c          write (iout,*) i,fdimless(i)
273         enddo
274 #endif
275       WRITE (iout,*) "Wchodze do call mysort1"
276         call mysort1(ncon,Fdimless,list_conf)
277       WRITE (iout,*) "Wychodze z call mysort1"
278 C#ifdef DEBUG
279         write (iout,*) "The FDIMLESS array after sorting"
280         do i=1,ncon
281           write (iout,*) i,list_conf(i),fdimless(i)
282         enddo
283 c#endif
284       WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
285         do i=1,ncon
286           totfree(i)=fdimless(i)
287         enddo
288         qfree=0.0d0
289         do i=1,ncon
290           qfree=qfree+exp(-fdimless(i)+fdimless(1))
291 c          write (iout,*) "fdimless", fdimless(i)
292         enddo
293 c        write (iout,*) "qfree",qfree
294         nlist=1
295         sumprob=0.0
296         write (iout,*) "ncon", ncon,maxstr_proc
297         do i=1,min0(ncon,maxstr_proc)-1 
298           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
299 C#ifdef DEBUG
300           write (iout,*) "tu szukaj ponizej 7"
301           write (iout,*) i,ib,beta_h(ib),
302      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
303      &     totfree(list_conf(i)),
304      &     -entfac(list_conf(i)),fdimless(i),sumprob
305 C#endif
306           if (sumprob.gt.prob_limit) goto 122
307 c          if (sumprob.gt.1.00d0) goto 122
308           nlist=nlist+1
309         enddo  
310   122   continue
311 #ifdef MPI
312       endif
313       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, 
314      &   IERROR)
315       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
316      &   IERROR)
317 c      do iproc=0,nprocs
318 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
319 c     &   " indend",indend(iproc) 
320 c      enddo
321       write (iout,*) "nlist",nlist
322 #endif
323       return
324       end
325 !--------------------------------------------------
326       subroutine mysort1(n, x, ipermut)
327       implicit none
328       integer i,j,imax,ipm,n
329       real x(n)
330       integer ipermut(n)
331       real xtemp
332       do i=1,n
333         xtemp=x(i)
334         imax=i
335         do j=i+1,n
336           if (x(j).lt.xtemp) then
337             imax=j
338             xtemp=x(j)
339           endif
340         enddo
341         x(imax)=x(i)
342         x(i)=xtemp
343         ipm=ipermut(imax)
344         ipermut(imax)=ipermut(i)
345         ipermut(i)=ipm
346       enddo
347       return
348       end