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