2 C Does modified entropic sampling in the space of minima.
3 implicit real*8 (a-h,o-z)
5 include 'COMMON.IOUNITS'
10 include 'COMMON.CHAIN'
13 include 'COMMON.CONTACTS'
14 include 'COMMON.CONTROL'
16 include 'COMMON.INTERACT'
17 include 'COMMON.THREAD'
18 include 'COMMON.NAMES'
19 logical accepted,not_done,over,ovrtim,error,lprint
20 integer MoveType,nbond
22 double precision RandOrPert
23 double precision varia(maxvar),elowest,ehighest,eold
24 double precision przes(3),obr(3,3)
25 double precision varold(maxvar)
27 double precision energia(0:n_ene),energia_ave(0:n_ene)
29 cd write (iout,*) 'print_mc=',print_mc
32 c---------------------------------------------------------------------------
33 C Initialize counters.
34 c---------------------------------------------------------------------------
35 C Total number of generated confs.
37 C Total number of moves. In general this won't be equal to the number of
38 C attempted moves, because we may want to reject some "bad" confs just by
41 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
46 C Initialize total and accepted number of moves of various kind.
51 C Total number of energy evaluations.
57 facee=1.0D0/(maxacc*delte)
59 C Read entropy from previous simulations.
61 read (ientin,*) indminn,indmaxx,emin,emax
62 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
65 entropy(i)=(emin+i*delte)*betbol
67 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
70 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
71 & ' emin=',emin,' emax=',emax
72 write (iout,'(/a)') 'Initial entropy'
74 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
77 C Read the pool of conformations
79 C----------------------------------------------------------------------------
80 C Entropy-sampling simulations with continually updated entropy
81 C Loop thru simulations
82 C----------------------------------------------------------------------------
84 C----------------------------------------------------------------------------
85 C Take a conformation from the pool
86 C----------------------------------------------------------------------------
92 write (iout,*) 'Took conformation',ii,' from the pool energy=',
94 call var_to_geom(nvar,varia)
95 C Print internal coordinates of the initial conformation
98 call gen_rand_conf(1,*20)
100 C----------------------------------------------------------------------------
101 C Compute and print initial energies.
102 C----------------------------------------------------------------------------
105 if (MyID.eq.MasterID) then
113 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
114 write (iout,'(/80(1h*)/a)') 'Initial energies:'
116 call etotal(energia(0))
118 call enerprint(energia(0))
119 C Minimize the energy of the first conformation.
121 call geom_to_var(nvar,varia)
122 call minimize(etot,varia,iretcode,nfun)
123 call etotal(energia(0))
125 write (iout,'(/80(1h*)/a/80(1h*))')
126 & 'Results of the first energy minimization:'
127 call enerprint(energia(0))
130 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
134 call contact(.false.,ncont,icont,co)
135 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
136 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
137 & 'RMS deviation from the reference structure:',rms,
138 & ' % of native contacts:',frac*100,' contact order:',co
139 write (istat,'(i5,11(1pe14.5))') 0,
140 & (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
142 write (istat,'(i5,9(1pe14.5))') 0,
143 & (energia(print_order(i)),i=1,nprint_ene),etot
146 neneval=neneval+nfun+1
147 if (.not. ent_read) then
148 C Initialize the entropy array
149 do i=-max_ene,max_ene
151 C Uncomment the line below for actual entropic sampling (start with uniform
152 C energy distribution).
154 C Uncomment the line below for multicanonical sampling (start with Boltzmann
156 entropy(i)=(emin+i*delte)*betbol
160 write (iout,'(/a)') 'Initial entropy'
162 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
166 call recv_stop_sig(Kwita)
167 if (whatsup.eq.1) then
168 call send_stop_sig(-2)
170 else if (whatsup.le.-2) then
172 else if (whatsup.eq.2) then
178 not_done = (iretcode.ne.11)
180 write (iout,'(/80(1h*)/20x,a/80(1h*))')
181 & 'Enter Monte Carlo procedure.'
183 call briefout(0,etot)
188 indeold=(eold-emin)/delte
189 deix=eold-(emin+indeold*delte)
190 dent=entropy(indeold+1)-entropy(indeold)
191 cd write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
192 cd write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
194 sold=entropy(indeold)+(dent/delte)*deix
196 write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
197 write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,
199 if (minim) call zapis(varia,etot)
201 C NACC is the counter for the accepted conformations of a given processor
203 C NACC_TOT counts the total number of accepted conformations
206 if (MyID.eq.MasterID) then
207 call receive_MCM_info
209 call send_MCM_info(2)
212 do iene=indminn,indmaxx
219 c----------------------------------------------------------------------------
225 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
226 & 'Beginning iteration #',it
227 C Initialize local counter.
228 ntrial=0 ! # of generated non-overlapping confs.
229 noverlap=0 ! # of overlapping confs.
231 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
233 C Retrieve the angles of previously accepted conformation
237 cd write (iout,'(a)') 'Old variables:'
238 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
239 call var_to_geom(nvar,varia)
245 C Decide whether to generate a random conformation or perturb the old one
246 RandOrPert=ran_number(0.0D0,1.0D0)
247 if (RandOrPert.gt.RanFract) then
249 & write (iout,'(a)') 'Perturbation-generated conformation.'
250 call perturb(error,lprint,MoveType,nbond,1.0D0)
252 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
253 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
254 & MoveType,' returned from PERTURB.'
261 nstart_grow=iran_num(3,nres)
263 & write (iout,'(2a,i3)') 'Random-generated conformation',
264 & ' - chain regrown from residue',nstart_grow
265 call gen_rand_conf(nstart_grow,*30)
267 call geom_to_var(nvar,varia)
268 cd write (iout,'(a)') 'New variables:'
269 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
271 if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)')
272 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
273 if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)')
274 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
275 call etotal(energia(0))
277 c call enerprint(energia(0))
278 c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
279 if (etot-elowest.gt.overlap_cut) then
280 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,
281 & ' Overlap detected in the current conf.; energy is',etot
285 if (noverlap.gt.maxoverlap) then
286 write (iout,'(a)') 'Too many overlapping confs.'
291 call minimize(etot,varia,iretcode,nfun)
292 cd write (iout,'(a)') 'Variables after minimization:'
293 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
294 call etotal(energia(0))
296 neneval=neneval+nfun+1
298 if (print_mc.gt.2) then
299 write (iout,'(a)') 'Total energies of trial conf:'
300 call enerprint(energia(0))
301 else if (print_mc.eq.1) then
302 write (iout,'(a,i6,a,1pe16.6)')
303 & 'Trial conformation:',ngen,' energy:',etot
305 C--------------------------------------------------------------------------
307 C--------------------------------------------------------------------------
310 & call accepting(etot,eold,scur,sold,varia,varold,
315 if (elowest.gt.etot) elowest=etot
316 if (ehighest.lt.etot) ehighest=etot
317 moves_acc(MoveType)=moves_acc(MoveType)+1
318 if (MoveType.eq.1) then
319 nbond_acc(nbond)=nbond_acc(nbond)+1
321 C Check against conformation repetitions.
322 irep=conf_comp(varia,etot)
323 #if defined(AIX) || defined(PGI) || defined(CRAY)
324 open (istat,file=statname,position='append')
326 open (istat,file=statname,access='append')
329 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
331 & przes,obr,non_conv)
333 call contact(.false.,ncont,icont,co)
334 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
336 & write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
337 & 'RMS deviation from the reference structure:',rms,
338 & ' % of native contacts:',frac*100,' contact order:',co
340 & write (istat,'(i5,11(1pe14.5))') it,
341 & (energia(print_order(i)),i=1,nprint_ene),etot,
343 elseif (print_stat) then
344 write (istat,'(i5,10(1pe14.5))') it,
345 & (energia(print_order(i)),i=1,nprint_ene),etot
349 & call statprint(nacc,nfun,iretcode,etot,elowest)
350 C Print internal coordinates.
351 if (print_int) call briefout(nacc,etot)
353 if (MyID.ne.MasterID) then
354 call recv_stop_sig(Kwita)
355 cd print *,'Processor:',MyID,' STOP=',Kwita
357 call send_MCM_info(2)
359 call send_MCM_info(1)
363 C Store the accepted conf. and its energy.
371 cd write (iout,*) 'Accepted conformation:'
372 cd write (iout,*) (rad2deg*varia(i),i=1,nphi)
373 if (minim) call zapis(varia,etot)
375 ener(i,nsave)=energia(i)
377 ener(n_ene+1,nsave)=etot
378 ener(n_ene+2,nsave)=frac
380 nminima(irep)=nminima(irep)+1.0D0
381 c print *,'irep=',irep,' nminima=',nminima(irep)
383 if (Kwita.eq.0) call recv_stop_sig(kwita)
388 if (MyID.eq.MasterID) then
389 call receive_MCM_info
390 if (nacc_tot.ge.maxacc) accepted=.true.
393 if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
394 C Take a conformation from the pool
399 write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
401 & 'Take conformation',ii,' from the pool energy=',epool(ii)
403 & write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
405 endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
409 if (MyID.eq.MasterID) then
410 call receive_MCM_info
412 if (Kwita.eq.0) call recv_stop_sig(kwita)
414 if (ovrtim()) WhatsUp=-1
415 cd write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
416 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
418 cd write (iout,*) 'not_done=',not_done
421 print *,'Processor',MyID,
422 & ' has received STOP signal =',Kwita,' in EntSamp.'
423 cd print *,'not_done=',not_done
424 if (Kwita.lt.-1) WhatsUp=Kwita
425 else if (nacc_tot.ge.maxacc) then
426 print *,'Processor',MyID,' calls send_stop_sig,',
427 & ' because a sufficient # of confs. have been collected.'
428 cd print *,'not_done=',not_done
429 call send_stop_sig(-1)
430 else if (WhatsUp.eq.-1) then
431 print *,'Processor',MyID,
432 & ' calls send_stop_sig because of timeout.'
433 cd print *,'not_done=',not_done
434 call send_stop_sig(-2)
439 C-----------------------------------------------------------------
440 C... Construct energy histogram & update entropy
441 C-----------------------------------------------------------------
445 write (iout,*) 'Processor',MyID,
446 & ' is broadcasting ERROR-STOP signal.'
447 write (*,*) 'Processor',MyID,
448 & ' is broadcasting ERROR-STOP signal.'
449 call send_stop_sig(-3)
453 if (MyID.eq.MasterID) then
454 c call receive_MCM_results
455 call receive_energies
458 if (esave(i).lt.elowest) elowest=esave(i)
459 if (esave(i).gt.ehighest) ehighest=esave(i)
461 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
462 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
463 & ' Highest energy',ehighest
464 if (isweep.eq.1 .and. .not.ent_read) then
467 write (iout,*) 'EMAX=',emax
469 indmaxx=(ehighest-emin)/delte
472 do i=-max_ene,max_ene
473 entropy(i)=(emin+i*delte)*betbol
477 indmin=(elowest-emin)/delte
478 indmax=(ehighest-emin)/delte
479 if (indmin.lt.indminn) indminn=indmin
480 if (indmax.gt.indmaxx) indmaxx=indmax
482 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
483 C Construct energy histogram
485 inde=(esave(i)-emin)/delte
486 nhist(inde)=nhist(inde)+nminima(i)
488 C Update entropy (density of states)
490 if (nhist(i).gt.0) then
491 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
495 Cd entropy(i)=1.0D+10
497 write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
498 & 'End of macroiteration',isweep
499 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
500 & ' Ehighest=',ehighest
501 write (iout,'(a)') 'Frequecies of minima'
503 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
505 write (iout,'(/a)') 'Energy histogram'
507 write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
509 write (iout,'(/a)') 'Entropy'
511 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
513 C-----------------------------------------------------------------
514 C... End of energy histogram construction
515 C-----------------------------------------------------------------
517 entropy(-max_ene-4)=dfloat(indminn)
518 entropy(-max_ene-3)=dfloat(indmaxx)
519 entropy(-max_ene-2)=emin
520 entropy(-max_ene-1)=emax
522 cd print *,entname,ientout
523 open (ientout,file=entname,status='unknown')
524 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
526 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
530 write (iout,'(a)') 'Frequecies of minima'
532 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
534 c call send_MCM_results
536 call receive_MCM_update
537 indminn=entropy(-max_ene-4)
538 indmaxx=entropy(-max_ene-3)
539 emin=entropy(-max_ene-2)
540 emax=entropy(-max_ene-1)
541 write (iout,*) 'Received from master:'
542 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
543 & ' emin=',emin,' emax=',emax
544 write (iout,'(/a)') 'Entropy'
546 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
549 if (WhatsUp.lt.-1) return
551 if (ovrtim() .or. WhatsUp.lt.0) return
554 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
555 call statprint(nacc,nfun,iretcode,etot,elowest)
557 & 'Statistics of multiple-bond motions. Total motions:'
558 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
559 write (iout,'(a)') 'Accepted motions:'
560 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
561 write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
562 write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
564 C---------------------------------------------------------------------------
566 C---------------------------------------------------------------------------
570 if (isweep.eq.nsweep .and. it.ge.maxacc)
571 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
574 c------------------------------------------------------------------------------
575 subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
576 implicit real*8 (a-h,o-z)
580 include 'COMMON.IOUNITS'
583 include 'COMMON.INFO'
586 double precision ecur,eold,xx,ran_number,bol
587 double precision x(maxvar),xold(maxvar)
588 double precision tole /1.0D-1/, tola /5.0D0/
590 C Check if the conformation is similar.
591 cd write (iout,*) 'Enter ACCEPTING'
592 cd write (iout,*) 'Old PHI angles:'
593 cd write (iout,*) (rad2deg*xold(i),i=1,nphi)
594 cd write (iout,*) 'Current angles'
595 cd write (iout,*) (rad2deg*x(i),i=1,nphi)
596 cd ddif=dif_ang(nphi,x,xold)
597 cd write (iout,*) 'Angle norm:',ddif
598 cd write (iout,*) 'ecur=',ecur,' emax=',emax
599 if (ecur.gt.emax) then
602 & write (iout,'(a)') 'Conformation rejected as too high in energy'
604 else if (dabs(ecur-eold).lt.tole .and.
605 & dif_ang(nphi,x,xold).lt.tola) then
608 & write (iout,'(a)') 'Conformation rejected as too similar'
611 C Else evaluate the entropy of the conf and compare it with that of the previous
613 indecur=(ecur-emin)/delte
614 if (iabs(indecur).gt.max_ene) then
615 write (iout,'(a,2i5)')
616 & 'Accepting: Index out of range:',indecur
618 else if (indecur.eq.indmaxx) then
619 scur=entropy(indecur)
620 if (print_mc.gt.0) write (iout,*)'Energy boundary reached',
621 & indmaxx,indecur,entropy(indecur)
623 deix=ecur-(emin+indecur*delte)
624 dent=entropy(indecur+1)-entropy(indecur)
625 scur=entropy(indecur)+(dent/delte)*deix
627 cd print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
628 cd & ' scur=',scur,' eold=',eold,' sold=',sold
629 cd print *,'deix=',deix,' dent=',dent,' delte=',delte
630 if (print_mc.gt.1) then
631 write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
632 write(iout,*)'eold=',eold,' sold=',sold
634 if (scur.le.sold) then
637 C Else carry out acceptance test
638 xx=ran_number(0.0D0,1.0D0)
640 if (xxh.gt.50.0D0) then
647 if (print_mc.gt.0) write (iout,'(a)')
648 & 'Conformation accepted.'
651 if (print_mc.gt.0) write (iout,'(a)')
652 & 'Conformation rejected.'
657 c-----------------------------------------------------------------------------
659 implicit real*8 (a-h,o-z)
661 include 'COMMON.IOUNITS'
666 double precision varia(maxvar)
667 print '(a)','Call READ_POOL'
670 read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
671 if (epool(npool).eq.0.0D0) goto 10
672 call read_angles(intin,*10)
673 call geom_to_var(nvar,xpool(1,npool))
677 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
678 if (print_mc.gt.2) then
680 write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',
682 write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
684 endif ! (print_mc.gt.2)