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