saxs and adam's corrections to multichain
[unres.git] / source / wham / src / make_ensemble1.F
1       subroutine make_ensembles(islice,*)
2 ! construct the conformational ensembles at REMD temperatures
3       implicit none
4       include "DIMENSIONS"
5       include "DIMENSIONS.ZSCOPT"
6       include "DIMENSIONS.FREE"
7 #ifdef MPI
8       include "mpif.h"
9       include "COMMON.MPI"
10       integer ierror,errcode,status(MPI_STATUS_SIZE) 
11 #endif
12       include "COMMON.IOUNITS"
13       include "COMMON.CONTROL"
14       include "COMMON.FREE"
15       include "COMMON.ENERGIES"
16       include "COMMON.FFIELD"
17       include "COMMON.INTERACT"
18       include "COMMON.SBRIDGE"
19       include "COMMON.CHAIN"
20       include "COMMON.PROTFILES"
21       include "COMMON.PROT"
22       real*4 csingle(3,maxres2)
23       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
24      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
25       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
26      &      escloc,ehomology_constr,
27      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
28      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt,
29      &      edfadis,edfator,edfanei,edfabet
30       integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
31       double precision qfree,sumprob,eini,efree,rmsdev
32       character*80 bxname
33       character*2 licz1,licz2
34       character*3 licz3,licz4
35       character*5 ctemper
36       integer ilen
37       external ilen
38       real*4 Fdimless(MaxStr),Fdimless_(MaxStr)
39       double precision enepot(MaxStr)
40       integer iperm(MaxStr)
41       integer islice
42  
43 #ifdef MPI
44       if (me.eq.Master) then
45 #endif
46       write (licz2,'(bz,i2.2)') islice
47       if (nslice.eq.1) then
48         if (.not.separate_parset) then
49           bxname = prefix(:ilen(prefix))//".bx"
50         else
51           write (licz3,'(bz,i3.3)') myparm
52           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
53         endif
54       else
55         if (.not.separate_parset) then
56           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
57         else
58           write (licz3,'(bz,i3.3)') myparm
59           bxname = prefix(:ilen(prefix))//"par_"//licz3//
60      &      "_slice_"//licz2//".bx"
61         endif
62       endif
63       open (ientout,file=bxname,status="unknown",
64      &  form="unformatted",access="direct",recl=lenrec1)
65 #ifdef MPI
66       endif
67 #endif
68       do iparm=1,nParmSet
69         if (iparm.ne.iparmprint) exit
70         call restore_parm(iparm)
71         do ib=1,nT_h(iparm)
72 #ifdef DEBUG
73           write (iout,*) "iparm",iparm," ib",ib
74 #endif
75           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
76 c          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
77 c          quotl=1.0d0
78 c          kfacl=1.0d0
79 c          do l=1,5
80 c            quotl1=quotl
81 c            quotl=quotl*quot
82 c            kfacl=kfacl*kfac
83 c            fT(l)=kfacl/(kfacl-1.0d0+quotl)
84 c          enddo
85             if (rescale_mode.eq.1) then
86               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
87 #if defined(FUNCTH)
88               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
89               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
90 #elif defined(FUNCT)
91               ft(6)=quot
92 #else
93               ft(6)=1.0d0
94 #endif
95               quotl=1.0d0
96               kfacl=1.0d0
97               do l=1,5
98                 quotl1=quotl
99                 quotl=quotl*quot
100                 kfacl=kfacl*kfac
101                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
102               enddo
103             else if (rescale_mode.eq.2) then
104               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
105 #if defined(FUNCTH)
106               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
107               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
108 #elif defined(FUNCT)
109               ft(6)=quot
110 #else 
111               ft(6)=1.0d0
112 #endif
113               quotl=1.0d0
114               do l=1,5
115                 quotl=quotl*quot
116                 fT(l)=1.12692801104297249644d0/
117      &             dlog(dexp(quotl)+dexp(-quotl))
118               enddo
119 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
120             else if (rescale_mode.eq.0) then
121               do l=1,5
122                 fT(l)=0.0d0
123               enddo
124             else
125               write (iout,*)
126      &        "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
127               call flush(iout)
128               return1
129             endif
130 #ifdef MPI
131           do i=1,scount(me1)
132 #else
133           do i=1,ntot(islice)
134 #endif
135             evdw=enetb(1,i,iparm)
136             evdw_t=enetb(21,i,iparm)
137 #ifdef SCP14
138             evdw2_14=enetb(17,i,iparm)
139             evdw2=enetb(2,i,iparm)+evdw2_14
140 #else
141             evdw2=enetb(2,i,iparm)
142             evdw2_14=0.0d0
143 #endif
144 #ifdef SPLITELE
145             ees=enetb(3,i,iparm)
146             evdw1=enetb(16,i,iparm)
147 #else
148             ees=enetb(3,i,iparm)
149             evdw1=0.0d0
150 #endif
151             ecorr=enetb(4,i,iparm)
152             ecorr5=enetb(5,i,iparm)
153             ecorr6=enetb(6,i,iparm)
154             eel_loc=enetb(7,i,iparm)
155             eello_turn3=enetb(8,i,iparm)
156             eello_turn4=enetb(9,i,iparm)
157             eturn6=enetb(10,i,iparm)
158             ebe=enetb(11,i,iparm)
159             escloc=enetb(12,i,iparm)
160             etors=enetb(13,i,iparm)
161             etors_d=enetb(14,i,iparm)
162             ehpb=enetb(15,i,iparm)
163             estr=enetb(18,i,iparm)
164             esccor=enetb(19,i,iparm)
165             edihcnstr=enetb(20,i,iparm)
166             ehomology_constr=enetb(22,i,iparm)
167             edfadis=enetb(23,i,iparm)
168             edfator=enetb(24,i,iparm)
169             edfanei=enetb(25,i,iparm)
170             edfabet=enetb(26,i,iparm)
171             if (homol_nset.gt.1)
172      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
173 #ifdef SPLITELE
174             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
175      &      +wvdwpp*evdw1
176      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
177      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
178      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
179      &      +ft(2)*wturn3*eello_turn3
180      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
181      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
182      &      +wbond*estr+ehomology_constr
183      &      +wdfa_dist*edfadis
184      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
185 #else
186             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
187      &      +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
194      &      +wbond*estr+ehomology_constr
195      &      +wdfa_dist*edfadis
196      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
197 #endif
198 #ifdef MPI
199             Fdimless_(i)=
200      &        beta_h(ib,iparm)*etot-entfac(i)
201             potE(i,iparm)=etot
202 #ifdef DEBUG
203             write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
204      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
205      &       -entfac(i),Fdimless(i)
206 #endif
207 #else
208             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
209             potE(i,iparm)=etot
210 #endif
211           enddo   ! i
212 #ifdef MPI
213           call MPI_Gatherv(Fdimless_(1),scount(me),
214      &     MPI_REAL,Fdimless(1),
215      &     scount(0),idispl(0),MPI_REAL,Master,
216      &     WHAM_COMM, IERROR)
217 #ifdef DEBUG
218           call MPI_Gatherv(potE(1,iparm),scount(me),
219      &     MPI_DOUBLE_PRECISION,potE(1,iparm),
220      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
221      &     WHAM_COMM, IERROR)
222           call MPI_Gatherv(entfac(1),scount(me),
223      &     MPI_DOUBLE_PRECISION,entfac(1),
224      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
225      &     WHAM_COMM, IERROR)
226 #endif
227           if (me.eq.Master) then
228 #ifdef DEBUG
229           write (iout,*) "The FDIMLESS array before sorting"
230           do i=1,ntot(islice)
231             write (iout,*) i,fdimless(i)
232           enddo
233 #endif
234 #endif
235           do i=1,ntot(islice)
236             iperm(i)=i
237           enddo
238           call mysort1(ntot(islice),Fdimless,iperm)
239 #ifdef DEBUG
240           write (iout,*) "The FDIMLESS array after sorting"
241           do i=1,ntot(islice)
242             write (iout,*) i,iperm(i),fdimless(i)
243           enddo
244 #endif
245           qfree=0.0d0
246           do i=1,ntot(islice)
247             qfree=qfree+exp(-fdimless(i)+fdimless(1))
248           enddo
249 c          write (iout,*) "qfree",qfree
250           nlist=1
251           sumprob=0.0
252           do i=1,min0(ntot(islice),ensembles) 
253             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
254 #ifdef DEBUG
255             write (iout,*) i,ib,beta_h(ib,iparm),
256      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
257      &       potE(iperm(i),iparm),
258      &       -entfac(iperm(i)),fdimless(i),sumprob
259 #endif
260             if (sumprob.gt.0.99d0) goto 122
261             nlist=nlist+1
262           enddo  
263   122     continue
264 #ifdef MPI
265           endif
266           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM, 
267      &       IERROR)
268           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
269      &       IERROR)
270           do i=1,nlist
271             ii=iperm(i)
272             iproc=0
273             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
274               iproc=iproc+1
275             enddo
276             if (iproc.ge.nprocs) then
277               write (iout,*) "Fatal error: processor out of range",iproc
278               call flush(iout)
279               if (bxfile) then
280                 close (ientout)
281               else
282                 close (ientout,status="delete")
283               endif
284               return1
285             endif
286             ik=ii-indstart(iproc)+1
287             if (iproc.ne.Master) then
288               if (me.eq.iproc) then
289 #ifdef DEBUG
290                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
291      &           " energy",potE(ik,iparm)
292 #endif
293                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
294      &            Master,i,WHAM_COMM,IERROR)
295               else if (me.eq.Master) then
296                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
297      &            WHAM_COMM,STATUS,IERROR)
298               endif
299             else if (me.eq.Master) then
300               enepot(i)=potE(ik,iparm)
301             endif
302           enddo
303 #else
304           do i=1,nlist
305             enepot(i)=potE(iperm(i),iparm)
306           enddo
307 #endif
308 #ifdef MPI
309           if (me.eq.Master) then
310 #endif
311           write(licz3,'(bz,i3.3)') iparm
312           write(licz2,'(bz,i2.2)') islice
313           if (temper.lt.100.0d0) then
314             write(ctemper,'(f3.0)') temper
315           else if (temper.lt.1000.0) then
316             write (ctemper,'(f4.0)') temper
317           else
318             write (ctemper,'(f5.0)') temper
319           endif
320           if (nparmset.eq.1) then
321             if (separate_parset) then
322               write(licz4,'(bz,i3.3)') myparm
323               pdbname=prefix(:ilen(prefix))//"_par"//licz4
324             else
325               pdbname=prefix(:ilen(prefix))
326             endif
327           else
328             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
329           endif
330           if (nslice.eq.1) then
331             pdbname=pdbname(:ilen(pdbname))//"_T_"//
332      &        ctemper(:ilen(ctemper))//"pdb"
333           else
334             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
335      &        ctemper(:ilen(ctemper))//"pdb" 
336           endif
337           open(ipdb,file=pdbname)
338           do i=1,nlist
339             read (ientout,rec=iperm(i))
340      &        ((csingle(l,k),l=1,3),k=1,nres),
341      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
342      &        nss,(ihpb(k),jhpb(k),k=1,nss),
343      &        eini,efree,rmsdev,iscor
344             do j=1,2*nres
345               do k=1,3
346                 c(k,j)=csingle(k,j)
347               enddo
348             enddo
349             eini=fdimless(i)
350             call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
351 cc          if (temper.eq.300.0d0) then
352 cc          write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree
353 cc          flush(iout)
354 cc          endif
355           enddo
356 cc          if (temper.eq.300.0d0) then
357 cc          write(iout,*) 'Gral i=',i ,eini,enepot(i),efree
358 cc          flush(iout)
359 cc          endif
360 #ifdef MPI
361           endif
362 #endif
363         enddo     ! ib
364       enddo       ! iparm
365       if (bxfile) then
366         close(ientout)
367       else
368         close(ientout,status="delete")
369       endif
370       return
371       end
372 !--------------------------------------------------
373       subroutine mysort1(n, x, ipermut)
374       implicit none
375       integer i,j,imax,ipm,n
376       real x(n)
377       integer ipermut(n)
378       real xtemp
379       do i=1,n
380         xtemp=x(i)
381         imax=i
382         do j=i+1,n
383           if (x(j).lt.xtemp) then
384             imax=j
385             xtemp=x(j)
386           endif
387         enddo
388         x(imax)=x(i)
389         x(i)=xtemp
390         ipm=ipermut(imax)
391         ipermut(imax)=ipermut(i)
392         ipermut(i)=ipm
393       enddo
394       return
395       end