Fixed the "***" issue in PDB output format by shifting the center of each molecule...
[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       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             ehomology_constr=enetb(22,i,iparm)
166             if (homol_nset.gt.1)
167      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
168 #ifdef SPLITELE
169             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
170      &      +wvdwpp*evdw1
171      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
172      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
173      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
174      &      +ft(2)*wturn3*eello_turn3
175      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
176      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
177      &      +wbond*estr+ehomology_constr
178 #else
179             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
180      &      +ft(1)*welec*(ees+evdw1)
181      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
182      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
183      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
184      &      +ft(2)*wturn3*eello_turn3
185      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
186      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
187      &      +wbond*estr+ehomology_constr
188 #endif
189 #ifdef MPI
190             Fdimless(i)=
191      &        beta_h(ib,iparm)*etot-entfac(i)
192             potE(i,iparm)=etot
193 #ifdef DEBUG
194             write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
195      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
196      &       -entfac(i),Fdimless(i)
197 #endif
198 #else
199             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
200             potE(i,iparm)=etot
201 #endif
202           enddo   ! i
203 #ifdef MPI
204           call MPI_Gatherv(Fdimless(1),scount(me),
205      &     MPI_REAL,Fdimless(1),
206      &     scount(0),idispl(0),MPI_REAL,Master,
207      &     WHAM_COMM, IERROR)
208 #ifdef DEBUG
209           call MPI_Gatherv(potE(1,iparm),scount(me),
210      &     MPI_DOUBLE_PRECISION,potE(1,iparm),
211      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
212      &     WHAM_COMM, IERROR)
213           call MPI_Gatherv(entfac(1),scount(me),
214      &     MPI_DOUBLE_PRECISION,entfac(1),
215      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
216      &     WHAM_COMM, IERROR)
217 #endif
218           if (me.eq.Master) then
219 #ifdef DEBUG
220           write (iout,*) "The FDIMLESS array before sorting"
221           do i=1,ntot(islice)
222             write (iout,*) i,fdimless(i)
223           enddo
224 #endif
225 #endif
226           do i=1,ntot(islice)
227             iperm(i)=i
228           enddo
229           call mysort1(ntot(islice),Fdimless,iperm)
230 #ifdef DEBUG
231           write (iout,*) "The FDIMLESS array after sorting"
232           do i=1,ntot(islice)
233             write (iout,*) i,iperm(i),fdimless(i)
234           enddo
235 #endif
236           qfree=0.0d0
237           do i=1,ntot(islice)
238             qfree=qfree+exp(-fdimless(i)+fdimless(1))
239           enddo
240 c          write (iout,*) "qfree",qfree
241           nlist=1
242           sumprob=0.0
243           do i=1,min0(ntot(islice),ensembles) 
244             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
245 #ifdef DEBUG
246             write (iout,*) i,ib,beta_h(ib,iparm),
247      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
248      &       potE(iperm(i),iparm),
249      &       -entfac(iperm(i)),fdimless(i),sumprob
250 #endif
251             if (sumprob.gt.0.99d0) goto 122
252             nlist=nlist+1
253           enddo  
254   122     continue
255 #ifdef MPI
256           endif
257           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM, 
258      &       IERROR)
259           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
260      &       IERROR)
261           do i=1,nlist
262             ii=iperm(i)
263             iproc=0
264             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
265               iproc=iproc+1
266             enddo
267             if (iproc.ge.nprocs) then
268               write (iout,*) "Fatal error: processor out of range",iproc
269               call flush(iout)
270               if (bxfile) then
271                 close (ientout)
272               else
273                 close (ientout,status="delete")
274               endif
275               return1
276             endif
277             ik=ii-indstart(iproc)+1
278             if (iproc.ne.Master) then
279               if (me.eq.iproc) then
280 #ifdef DEBUG
281                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
282      &           " energy",potE(ik,iparm)
283 #endif
284                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
285      &            Master,i,WHAM_COMM,IERROR)
286               else if (me.eq.Master) then
287                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
288      &            WHAM_COMM,STATUS,IERROR)
289               endif
290             else if (me.eq.Master) then
291               enepot(i)=potE(ik,iparm)
292             endif
293           enddo
294 #else
295           do i=1,nlist
296             enepot(i)=potE(iperm(i),iparm)
297           enddo
298 #endif
299 #ifdef MPI
300           if (me.eq.Master) then
301 #endif
302           write(licz3,'(bz,i3.3)') iparm
303           write(licz2,'(bz,i2.2)') islice
304           if (temper.lt.100.0d0) then
305             write(ctemper,'(f3.0)') temper
306           else if (temper.lt.1000.0) then
307             write (ctemper,'(f4.0)') temper
308           else
309             write (ctemper,'(f5.0)') temper
310           endif
311           if (nparmset.eq.1) then
312             if (separate_parset) then
313               write(licz4,'(bz,i3.3)') myparm
314               pdbname=prefix(:ilen(prefix))//"_par"//licz4
315             else
316               pdbname=prefix(:ilen(prefix))
317             endif
318           else
319             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
320           endif
321           if (nslice.eq.1) then
322             pdbname=pdbname(:ilen(pdbname))//"_T_"//
323      &        ctemper(:ilen(ctemper))//"pdb"
324           else
325             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
326      &        ctemper(:ilen(ctemper))//"pdb" 
327           endif
328           open(ipdb,file=pdbname)
329           do i=1,nlist
330             read (ientout,rec=iperm(i))
331      &        ((csingle(l,k),l=1,3),k=1,nres),
332      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
333      &        nss,(ihpb(k),jhpb(k),k=1,nss),
334      &        eini,efree,rmsdev,iscor
335             do j=1,2*nres
336               do k=1,3
337                 c(k,j)=csingle(k,j)
338               enddo
339             enddo
340             eini=fdimless(i)
341             call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
342 cc          if (temper.eq.300.0d0) then
343 cc          write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree
344 cc          flush(iout)
345 cc          endif
346           enddo
347 cc          if (temper.eq.300.0d0) then
348 cc          write(iout,*) 'Gral i=',i ,eini,enepot(i),efree
349 cc          flush(iout)
350 cc          endif
351 #ifdef MPI
352           endif
353 #endif
354         enddo     ! ib
355       enddo       ! iparm
356       if (bxfile) then
357         close(ientout)
358       else
359         close(ientout,status="delete")
360       endif
361       return
362       end
363 !--------------------------------------------------
364       subroutine mysort1(n, x, ipermut)
365       implicit none
366       integer i,j,imax,ipm,n
367       real x(n)
368       integer ipermut(n)
369       real xtemp
370       do i=1,n
371         xtemp=x(i)
372         imax=i
373         do j=i+1,n
374           if (x(j).lt.xtemp) then
375             imax=j
376             xtemp=x(j)
377           endif
378         enddo
379         x(imax)=x(i)
380         x(i)=xtemp
381         ipm=ipermut(imax)
382         ipermut(imax)=ipermut(i)
383         ipermut(i)=ipm
384       enddo
385       return
386       end