1 subroutine proc_data(nvarr,x,*)
4 include "DIMENSIONS.ZSCOPT"
7 integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
10 include "COMMON.CONTROL"
11 include "COMMON.ENERGIES"
12 include "COMMON.VMCPAR"
13 include "COMMON.OPTIM"
14 include "COMMON.WEIGHTS"
15 include "COMMON.PROTNAME"
16 include "COMMON.IOUNITS"
17 include "COMMON.FFIELD"
18 include "COMMON.PROTFILES"
19 include "COMMON.COMPAR"
20 include "COMMON.CHAIN"
21 include "COMMON.CLASSES"
22 include "COMMON.THERMAL"
23 include "COMMON.WEIGHTDER"
25 include "COMMON.SBRIDGE"
26 include "COMMON.INTERACT"
27 double precision fac0,t0
28 double precision fT(5),ftprim(5),ftbis(5),quot,quotl,quotl1,
29 & denom,kfac,kfacl,logfac,eplus,eminus,tanhT,entfacmax,
31 real csingle(3,maxres2)
32 double precision creff(3,maxres2),cc(3,maxres2)
33 double precision przes(3),obrot(3,3),etoti
36 integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng,
37 & inn,in,iref,ib,ibatch,num
38 integer l1,l2,ncon_all
39 double precision rcyfr
40 integer jj,istart_conf,iend_conf,ipass_conf,iii
43 double precision ej,emaxj,sumP,sumP_t
45 double precision x(max_paropt)
46 double precision wconf(maxstr_proc)
49 double precision rms,rmsnat,qwolynes,qwolyness
53 write (iout,*) "fac0",fac0
55 c Determine the number of variables and the initial vector of optimizable
59 if (restart_flag) then
60 call w2x(nvarr,x_orig,*11)
64 open (88,file=prefix(:ilen(prefix))//'.restin',status='unknown')
65 read(88,*,end=99,err=99) (x(i),i=1,nvarr)
69 & "Variables replaced with values from the restart file."
71 write (iout,'(i5,f15.5)') i,x(i)
77 & "Error - restart file nonexistent, incompatible or damaged."
79 call MPI_Finalize(Ierror)
82 & "Error - restart file nonexistent, incompatible or damaged."
94 c Setup weights for UNRES
123 call restore_molinfo(iprot)
125 DO IBATCH=1,natlike(iprot)+2
127 c Calculate temperature-dependent weight scaling factors
128 do ib=1,nbeta(ibatch,iprot)
129 fac=betaT(ib,ibatch,iprot)
131 if (rescale_mode.eq.0) then
137 else if (rescale_mode.eq.1) then
144 denom=kfacl-1.0d0+quotl
146 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
147 ftbis(l)=l*kfacl*quotl1*
148 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
149 write (iout,*) "ib",ib," beta",betaT(ib,ibatch,iprot),
150 & " quot",quot," l",l,fT(l),fTprim(l),fTbis(l)
152 else if (rescale_mode.eq.2) then
159 logfac=1.0d0/dlog(eplus+eminus)
160 tanhT=(eplus-eminus)/(eplus+eminus)
161 fT(l)=1.12692801104297249644d0*logfac
162 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
163 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
164 & 2*l*quotl1/T0*logfac*
165 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
169 write (iout,*) "Unknown RESCALE_MODE",rescale_mode
174 escal(k,ib,ibatch,iprot)=1.0d0
175 escalprim(k,ib,ibatch,iprot)=0.0d0
176 escalbis(k,ib,ibatch,iprot)=0.0d0
178 escal(3,ib,ibatch,iprot)=ft(1)
179 escal(4,ib,ibatch,iprot)=ft(3)
180 escal(5,ib,ibatch,iprot)=ft(4)
181 escal(6,ib,ibatch,iprot)=ft(5)
182 escal(7,ib,ibatch,iprot)=ft(2)
183 escal(8,ib,ibatch,iprot)=ft(2)
184 escal(9,ib,ibatch,iprot)=ft(3)
185 escal(10,ib,ibatch,iprot)=ft(5)
186 escal(13,ib,ibatch,iprot)=ft(1)
187 escal(14,ib,ibatch,iprot)=ft(2)
188 escal(19,ib,ibatch,iprot)=ft(1)
189 escalprim(3,ib,ibatch,iprot)=ftprim(1)
190 escalprim(4,ib,ibatch,iprot)=ftprim(3)
191 escalprim(5,ib,ibatch,iprot)=ftprim(4)
192 escalprim(6,ib,ibatch,iprot)=ftprim(5)
193 escalprim(7,ib,ibatch,iprot)=ftprim(2)
194 escalprim(8,ib,ibatch,iprot)=ftprim(2)
195 escalprim(9,ib,ibatch,iprot)=ftprim(3)
196 escalprim(10,ib,ibatch,iprot)=ftprim(5)
197 escalprim(13,ib,ibatch,iprot)=ftprim(1)
198 escalprim(14,ib,ibatch,iprot)=ftprim(2)
199 escalprim(19,ib,ibatch,iprot)=ftprim(1)
200 escalbis(3,ib,ibatch,iprot)=ftbis(1)
201 escalbis(4,ib,ibatch,iprot)=ftbis(3)
202 escalbis(5,ib,ibatch,iprot)=ftbis(4)
203 escalbis(6,ib,ibatch,iprot)=ftbis(5)
204 escalbis(7,ib,ibatch,iprot)=ftbis(2)
205 escalbis(8,ib,ibatch,iprot)=ftbis(2)
206 escalbis(9,ib,ibatch,iprot)=ftbis(3)
207 escalbis(10,ib,ibatch,iprot)=ftbis(5)
208 escalbis(13,ib,ibatch,iprot)=ftbis(1)
209 escalbis(14,ib,ibatch,iprot)=ftbis(2)
210 escalbis(19,ib,ibatch,iprot)=ftbis(1)
213 betmin(ibatch,iprot)=betaT(1,ibatch,iprot)
214 betmax(ibatch,iprot)=betaT(1,ibatch,iprot)
215 do ib=2,nbeta(ibatch,iprot)
216 if (betaT(ib,ibatch,iprot).lt.betmin(ibatch,iprot))
217 & betmin(ibatch,iprot)=betaT(ib,ibatch,iprot)
218 if (betaT(ib,ibatch,iprot).gt.betmax(ibatch,iprot))
219 & betmax(ibatch,iprot)=betaT(ib,ibatch,iprot)
221 write (iout,'(2(a,i5,2x),2(a,f7.2,2x))')
222 & "Protein",iprot," batch",ibatch,
223 & " betmin",betmin(ibatch,iprot)," betmax",betmax(ibatch,iprot)
229 c Compute energy and entropy histograms to filter out conformations
230 c with potntially too low or too high weights.
232 call determine_histograms
234 c Make the working list of conformations
236 call make_list(.true.,.true.,nvarr,x)
237 c 3/8/2013 Calculate RMSDs/Qs and the Ptab array
240 call restore_molinfo(iprot)
242 open (icbase,file=bprotfiles(iprot),status="old",
243 & form="unformatted",access="direct",recl=lenrec(iprot))
244 entfacmax=entfac(1,iprot)
245 sument=entfac(1,iprot)
246 sumentsq=entfac(1,iprot)**2
247 do i=2,ntot_work(iprot)
248 if (entfac(i,iprot).gt.entfacmax) entfacmax=entfac(i,iprot)
249 sument=sument+entfac(i,iprot)
250 sumentsq=sumentsq+entfac(i,iprot)**2
252 sument=sument/ntot_work(iprot)
253 sumentsq=dsqrt(sumentsq/ntot_work(iprot)-sument**2)
254 write (2,*) "entfacmax",entfacmax," entave",sument,
258 nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
262 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
263 ipass_conf=ipass_conf+1
264 write (iout,*) "Pass",ipass_conf
266 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
268 nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
271 do i=1,ntot(iprot),maxstr_proc
272 ipass_conf=ipass_conf+1
273 write (iout,*) "Pass",ipass_conf
275 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
278 c Read the chunk of conformations off a DA scratchfile.
280 call daread_ccoords(iprot,istart_conf,iend_conf)
282 IF (GEOM_AND_WHAM_WEIGHTS) THEN
285 do ib=1,nbeta(1,iprot)
286 write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
288 do iii=istart_conf,iend_conf
294 etoti=etoti+ww(k)*escal(k,ib,1,iprot)
297 e_total_T(iii,ib)=entfac(kk,iprot)
298 & +betaT(ib,1,iprot)*(etoti-elowest(ib,1,iprot))
302 write (iout,*) "Energies before Gather"
303 do iii=1,ntot_work(iprot)
304 write (iout,'(i5,100E12.5)') iii,
305 & (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
309 do ib=1,nbeta(1,iprot)
310 call MPI_AllGatherv(e_total_T(indstart(me1,iprot),ib),
311 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
312 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
316 write (iout,*) "Energies after Gather"
317 do iii=1,ntot_work(iprot)
318 write (iout,'(i5,100E12.5)') iii,
319 & (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
324 ENDIF ! GEOM_AND_WHAM_WEIGHTS
328 do iii=istart_conf,iend_conf
332 call restore_ccoords(iprot,ii)
333 c Calculate the rmsds of the current conformations and convert them into weights
334 c to approximate equal weighting of all points of the conformational space
340 do ib=1,nbeta(1,iprot)
343 do iref=1,ntot_work(iprot)
344 if (iref.ne.iii) then
345 read(icbase,rec=iref) ((csingle(l,k),l=1,3),k=1,nres),
346 & ((csingle(l,k),l=1,3),k=nnt+nres,nct+nres),
347 & nss,(ihpb(k),jhpb(k),k=1,nss),eini(iref,iprot),
348 & entfac(iref,iprot),
349 & ((nu(k,l,ii,iprot),k=1,maxdimnat),l=1,natlike(iprot))
357 write (iout,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
361 if (rmscomp(iprot)) then
362 call fitsq(rms,c(1,nstart_sup(iprot)),
363 & creff(1,nstart_sup(iprot)),
364 & nsup(iprot),przes,obrot,non_conv)
366 print *,'Error: FITSQ non-convergent, iref',iref
368 else if (rms.lt.-1.0d-6) then
369 print *,'Error: rms^2 = ',rms,iref
371 else if (rms.lt.0) then
377 write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
378 & e_total_T(iref,ib),
379 & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
380 & " rms",rms," fac",dexp(-0.5d0*(rms/sigma2(iprot))**2),
381 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
382 & dexp(-0.5d0*(rms/sigma2(iprot))**2)*
383 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
385 aux=dexp(-0.5d0*(rms/sigma2(iprot))**2)
387 rms=qwolyness(creff,iprot)
389 write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
390 & e_total_T(iref,ib),
391 & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
392 & " rms",-dlog(rms)," fac",rms,
393 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
394 & rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
401 do ib=1,nbeta(1,iprot)
402 if (GEOM_AND_WHAM_WEIGHTS) then
403 weirms(iii,ib) = weirms(iii,ib)
404 & +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
406 weirms(iii,ib) = weirms(iii,ib)+aux
412 write (iout,*) "weirms"
413 do iii=istart_conf,iend_conf
414 write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1,
422 nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
423 DO IB=1,NBETA(1,IPROT)
427 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
428 ipass_conf=ipass_conf+1
429 write (iout,*) "Pass",ipass_conf
431 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
433 nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
436 do i=1,ntot(iprot),maxstr_proc
437 ipass_conf=ipass_conf+1
438 write (iout,*) "Pass",ipass_conf
440 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
444 c Read the chunk of conformations off a DA scratchfile.
446 call daread_ccoords(iprot,istart_conf,iend_conf)
447 write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
448 do iii=istart_conf,iend_conf
452 call restore_ccoords(iprot,ii)
453 Ptab(jj,ib,iprot)=0.0d0
454 if (rmscomp(iprot)) then
455 do iref=1,nref(ib,iprot)
457 write (2,*) "nstart_sup",nstart_sup(iprot),
458 & " nend_sup",nend_sup(iprot)
459 do k=nstart_sup(iprot),nend_sup(iprot)
460 write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
461 & (cref(l,k,iref,ib,iprot),l=1,3)
463 do k=nstart_sup(iprot),nend_sup(iprot)
464 write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3),
465 & (cref(l,k+nres,iref,ib,iprot),l=1,3)
468 rms=rmsnat(ib,jj,iref,iprot)
470 write (iout,*) "iii",iii," ii",ii," jj",jj," rms",rms,
471 & "efree",efreeref(iref,ib,iprot)
473 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)
474 & +dexp(-efreeref(iref,ib,iprot)
475 & -0.5d0*(rms/sigma2(iprot))**2)
478 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
479 c search of conformational space.
480 aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
481 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
483 write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
484 c & " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot)
485 & "entfac",entfac(k,iprot)," weirms",weirms(k,iprot),
486 & " Ptab",Ptab(jj,ib,iprot)
488 #elif defined(WEIDIST)
489 write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
490 & " weirms",weirms(iii,ib)," Ptab/weirms",
491 & Ptab(jj,ib,iprot)/weirms(iii,ib)
492 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
493 #elif defined(WEICLUST)
494 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot)
496 sumP=sumP+Ptab(jj,ib,iprot)
498 do iref=1,nref(ib,iprot)
499 rms = qwolynes(ib,iref,iprot)
501 write (iout,*) "iii",iii," ii",ii," jj",jj,
502 & "iref",iref," rms",rms,-dlog(rms),
503 & "efree",efreeref(iref,ib,iprot)
505 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms*
506 & *dexp(-efreeref(iref,ib,iprot))
509 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
510 c search of conformational space.
511 aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
512 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
514 write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
515 & " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot)
517 #elif defined(WEICLUST)
518 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot)
519 #elif defined(WEIDIST)
520 write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
521 & " weirms",weirms(iii,ib)," Ptab/weirms",
522 & Ptab(jj,ib,iprot)/weirms(iii,ib)
523 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
525 sumP=sumP+Ptab(jj,ib,iprot)
528 write (iout,*) "ii",ii," jj",jj," ib",ib,
529 & "Ptab",Ptab(jj,ib,iprot)," sumP",sumP
534 write (iout,*) "ib",ib," sumP before REDUCE",sumP
537 call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION,
538 & MPI_SUM,Comm1,IERROR)
542 write (iout,*) "ib",ib," sumP after REDUCE",sumP
545 do iii=indstart(me1,iprot),indend(me1,iprot)
548 write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
550 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP
552 write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
560 write (iout,*) "ELOWEST at the end of PROC_DATA"
562 do ibatch=1,natlike(iprot)+2
563 do ib=1,nbeta(ibatch,iprot)
564 write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
565 & " elowest",elowest(ib,ibatch,iprot)
574 c-------------------------------------------------------------------------------
575 subroutine determine_histograms
578 include "DIMENSIONS.ZSCOPT"
580 parameter (maxind=1000)
583 integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
586 include "COMMON.COMPAR"
587 include "COMMON.ENERGIES"
588 include "COMMON.VMCPAR"
589 include "COMMON.WEIGHTS"
590 include "COMMON.PROTNAME"
591 include "COMMON.IOUNITS"
592 include "COMMON.FFIELD"
593 include "COMMON.PROTFILES"
594 include "COMMON.CHAIN"
595 include "COMMON.CONTROL"
596 include "COMMON.CLASSES"
597 integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt
598 double precision eminh,emaxh,tminh,tmaxh
599 integer he(0:maxind),
600 & ht(0:maxind),hemax,htmax
602 c Construct energy and entropy histograms
604 write (iout,*) "Calling DETERMINE HISTOGRAMS"
617 if (eini(iii,iprot).lt.eminh)
618 & eminh=eini(iii,iprot)
619 if (eini(iii,iprot).gt.emaxh)
620 & emaxh=eini(iii,iprot)
621 if (entfac(iii,iprot).lt.tminh)
622 & tminh=entfac(iii,iprot)
623 if (entfac(iii,iprot).gt.tmaxh)
624 & tmaxh=entfac(iii,iprot)
626 c write (2,*) "DETERMINE HISTOGRAMS 1"
630 inde=eini(iii,iprot)-eminh
631 indt=entfac(iii,iprot)-tminh
632 if (inde.le.maxind) he(inde)=he(inde)+1
633 if (indt.le.maxind) ht(indt)=ht(indt)+1
635 c write (2,*) "DETERMINE HISTOGRAMS 2"
639 c write (iout,*) "idelte",idelte," ideltt",ideltt
642 c write (iout,*) "i",i," he",he(i)," hemax",hemax
644 if (he(i).gt.hemax) hemax=he(i)
648 c write (iout,*) "i",i," ht",ht(i)," htmax",htmax
650 if (ht(i).gt.htmax) htmax=ht(i)
652 c write (2,*) "DETERMINE HISTOGRAMS 3"
654 hemax=hemax/hefac(iprot)
655 if (hemax.lt.hemax_low(iprot))
656 & hemax=hemax_low(iprot)
657 htmax=htmax/htfac(iprot)
658 if (htmax.lt.htmax_low(iprot))
659 & htmax=htmax_low(iprot)
661 do while (i.lt.idelte .and. he(i).lt.hemax)
664 c write (2,*) "DETERMINE HISTOGRAMS 4"
666 e_lowb(iprot)=eminh+i
668 do while (i.lt.ideltt .and. ht(i).lt.htmax)
671 t_lowb(iprot)=tminh+i
673 write (iout,*) "protein",iprot
674 write (iout,*) "energy histogram"
676 write (iout,'(f15.5,i10)') eminh+i,he(i)
678 write (iout,*) "entropy histogram"
680 write (iout,'(f15.5,i10)') tminh+i,ht(i)
682 write (iout,*) "emin",eminh," tmin",tminh,
683 & " hemax",hemax," htmax",htmax,
684 & " e_lowb",e_lowb(iprot),
685 & " t_lowb",t_lowb(iprot)
691 c------------------------------------------------------------------------------
692 subroutine imysort(n, x, ipermut)
702 if (x(j).lt.xtemp) then
710 ipermut(imax)=ipermut(i)
715 c--------------------------------------------------------------------------
716 subroutine write_conf_count
719 include "DIMENSIONS.ZSCOPT"
720 include "COMMON.CLASSES"
721 include "COMMON.COMPAR"
722 include "COMMON.VMCPAR"
723 include "COMMON.PROTNAME"
724 include "COMMON.IOUNITS"
725 include "COMMON.PROTFILES"
726 integer iprot,ibatch,i,ii,j,kk
731 write (iout,*)"Numbers of conformations after applying the cutoff"
734 write (iout,'(a,2x,a,i10)') "Protein",
735 & protname(iprot)(:ilen(protname(iprot))),
740 c-------------------------------------------------------------------------
741 double precision function fdum()