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