64e8c505cabe40df80c76fb6a9e1a8eb16130d1e
[unres.git] / source / cluster / wham / src / 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(5),fTprim(5),fTbis(5),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       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
25       double precision qfree,sumprob,eini,efree,rmsdev
26       character*80 bxname
27       character*2 licz1
28       character*5 ctemper
29       integer ilen
30       external ilen
31       real*4 Fdimless(maxconf)
32       double precision energia(0:max_ene)
33       do i=1,ncon
34         list_conf(i)=i
35       enddo
36 c      do i=1,ncon
37 c        write (iout,*) i,list_conf(i)
38 c      enddo
39 #ifdef MPI
40       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
41       call daread_ccoords(indstart(me),indend(me))
42 #endif
43 c      write (iout,*) "ncon",ncon
44       temper=1.0d0/(beta_h(ib)*1.987D-3)
45 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
46 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
47 c      quotl=1.0d0
48 c      kfacl=1.0d0
49 c      do l=1,5
50 c        quotl1=quotl
51 c        quotl=quotl*quot
52 c        kfacl=kfacl*kfac
53 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
54 c      enddo
55             if (rescale_mode.eq.1) then
56               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
57               quotl=1.0d0
58               kfacl=1.0d0
59               do l=1,5
60                 quotl1=quotl
61                 quotl=quotl*quot
62                 kfacl=kfacl*kfac
63                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
64               enddo
65             else if (rescale_mode.eq.2) then
66               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
67               quotl=1.0d0
68               do l=1,5
69                 quotl=quotl*quot
70                 fT(l)=1.12692801104297249644d0/
71      &             dlog(dexp(quotl)+dexp(-quotl))
72               enddo
73 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
74               call flush(iout)
75             endif
76
77 #ifdef MPI
78       do i=1,scount(me)
79         ii=i+indstart(me)-1
80 #else
81       do i=1,ncon
82         ii=i
83 #endif
84 c        write (iout,*) "i",i," ii",ii
85         call flush(iout)
86         if (ib.eq.1) then
87           do j=1,nres
88             do k=1,3
89               c(k,j)=allcart(k,j,i)
90               c(k,j+nres)=allcart(k,j+nres,i)
91             enddo
92           enddo
93           do k=1,3
94             c(k,nres+1)=c(k,1)
95             c(k,nres+nres)=c(k,nres)
96           enddo
97           nss=nss_all(i)
98           do j=1,nss
99             ihpb(j)=ihpb_all(j,i)
100             jhpb(j)=jhpb_all(j,i)
101           enddo 
102           call int_from_cart1(.false.)
103           call etotal(energia(0),fT)
104           totfree(i)=energia(0)
105 #ifdef DEBUG
106           write (iout,*) i," energia",(energia(j),j=0,19)
107 #endif
108           do k=1,max_ene
109             enetb(k,i)=energia(k)
110           enddo
111         endif
112         evdw=enetb(1,i)
113 #ifdef SCP14
114         evdw2_14=enetb(17,i)
115         evdw2=enetb(2,i)+evdw2_14
116 #else
117         evdw2=enetb(2,i)
118         evdw2_14=0.0d0
119 #endif
120 #ifdef SPLITELE
121         ees=enetb(3,i)
122         evdw1=enetb(16,i)
123 #else
124         ees=enetb(3,i)
125         evdw1=0.0d0
126 #endif
127         ecorr=enetb(4,i)
128         ecorr5=enetb(5,i)
129         ecorr6=enetb(6,i)
130         eel_loc=enetb(7,i)
131         eello_turn3=enetb(8,i)
132         eello_turn4=enetb(9,i)
133         eturn6=enetb(10,i)
134         ebe=enetb(11,i)
135         escloc=enetb(12,i)
136         etors=enetb(13,i)
137         etors_d=enetb(14,i)
138         ehpb=enetb(15,i)
139         estr=enetb(18,i)
140         esccor=enetb(19,i)
141         edihcnstr=enetb(20,i)
142 #ifdef SPLITELE
143         etot=wsc*evdw+wscp*evdw2+ft(1)*welec*ees+wvdwpp*evdw1
144      &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
145      &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
146      &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
147      &  +ft(2)*wturn3*eello_turn3
148      &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
149      &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
150      &  +wbond*estr
151 #else
152         etot=wsc*evdw+wscp*evdw2+ft(1)*welec*(ees+evdw1)
153      &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
154      &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
155      &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
156      &  +ft(2)*wturn3*eello_turn3
157      &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
158      &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
159      &  +wbond*estr
160 #endif
161         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
162         totfree(i)=etot
163 #ifdef DEBUG
164         write (iout,*) i,ii,ib,
165      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
166      &   entfac(ii),Fdimless(i)
167 #endif
168       enddo   ! i
169 #ifdef MPI
170       call MPI_Gatherv(Fdimless(1),scount(me),
171      & MPI_REAL,Fdimless(1),
172      & scount(0),idispl(0),MPI_REAL,Master,
173      & MPI_COMM_WORLD, IERROR)
174       call MPI_Gatherv(totfree(1),scount(me),
175      & MPI_DOUBLE_PRECISION,totfree(1),
176      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
177      & MPI_COMM_WORLD, IERROR)
178       call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
179      & MPI_DOUBLE_PRECISION,entfac(1),
180      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
181      & MPI_COMM_WORLD, IERROR)
182       if (me.eq.Master) then
183 #endif
184 #ifdef DEBUG
185         write (iout,*) "The FDIMLESS array before sorting"
186         do i=1,ncon
187           write (iout,*) i,fdimless(i)
188         enddo
189 #endif
190         call mysort1(ncon,Fdimless,list_conf)
191 #ifdef DEBUG
192         write (iout,*) "The FDIMLESS array after sorting"
193         do i=1,ncon
194           write (iout,*) i,list_conf(i),fdimless(i)
195         enddo
196 #endif
197         do i=1,ncon
198           totfree(i)=fdimless(i)
199         enddo
200         qfree=0.0d0
201         do i=1,ncon
202           qfree=qfree+exp(-fdimless(i)+fdimless(1))
203         enddo
204 c        write (iout,*) "qfree",qfree
205         nlist=1
206         sumprob=0.0
207         do i=1,min0(ncon,maxstr_proc)-1 
208           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
209 #ifdef DEBUG
210           write (iout,*) i,ib,beta_h(ib),
211      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
212      &     totfree(list_conf(i)),
213      &     -entfac(list_conf(i)),fdimless(i),sumprob
214 #endif
215           if (sumprob.gt.prob_limit) goto 122
216 c          if (sumprob.gt.1.00d0) goto 122
217           nlist=nlist+1
218         enddo  
219   122   continue
220 #ifdef MPI
221       endif
222       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, 
223      &   IERROR)
224       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
225      &   IERROR)
226 c      do iproc=0,nprocs
227 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
228 c     &   " indend",indend(iproc) 
229 c      enddo
230       write (iout,*) "nlist",nlist
231 #endif
232       return
233       end
234 !--------------------------------------------------
235       subroutine mysort1(n, x, ipermut)
236       implicit none
237       integer i,j,imax,ipm,n
238       real x(n)
239       integer ipermut(n)
240       real xtemp
241       do i=1,n
242         xtemp=x(i)
243         imax=i
244         do j=i+1,n
245           if (x(j).lt.xtemp) then
246             imax=j
247             xtemp=x(j)
248           endif
249         enddo
250         x(imax)=x(i)
251         x(i)=xtemp
252         ipm=ipermut(imax)
253         ipermut(imax)=ipermut(i)
254         ipermut(i)=ipm
255       enddo
256       return
257       end