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