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
100 call gen_rand_conf(jeden,nrestmp,*20)
102 C----------------------------------------------------------------------------
103 C Compute and print initial energies.
104 C----------------------------------------------------------------------------
107 if (MyID.eq.MasterID) then
115 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
116 write (iout,'(/80(1h*)/a)') 'Initial energies:'
118 call etotal(energia(0))
120 call enerprint(energia(0))
121 C Minimize the energy of the first conformation.
123 call geom_to_var(nvar,varia)
124 call minimize(etot,varia,iretcode,nfun)
125 call etotal(energia(0))
127 write (iout,'(/80(1h*)/a/80(1h*))')
128 & 'Results of the first energy minimization:'
129 call enerprint(energia(0))
133 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
137 call contact(.false.,ncont,icont,co)
138 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
139 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
140 & 'RMS deviation from the reference structure:',rms,
141 & ' % of native contacts:',frac*100,' contact order:',co
142 write (istat,'(i5,11(1pe14.5))') 0,
143 & (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
145 write (istat,'(i5,9(1pe14.5))') 0,
146 & (energia(print_order(i)),i=1,nprint_ene),etot
149 neneval=neneval+nfun+1
150 if (.not. ent_read) then
151 C Initialize the entropy array
152 do i=-max_ene,max_ene
154 C Uncomment the line below for actual entropic sampling (start with uniform
155 C energy distribution).
157 C Uncomment the line below for multicanonical sampling (start with Boltzmann
159 entropy(i)=(emin+i*delte)*betbol
163 write (iout,'(/a)') 'Initial entropy'
165 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
169 call recv_stop_sig(Kwita)
170 if (whatsup.eq.1) then
171 call send_stop_sig(-2)
173 else if (whatsup.le.-2) then
175 else if (whatsup.eq.2) then
181 not_done = (iretcode.ne.11)
183 write (iout,'(/80(1h*)/20x,a/80(1h*))')
184 & 'Enter Monte Carlo procedure.'
186 call briefout(0,etot)
191 indeold=(eold-emin)/delte
192 deix=eold-(emin+indeold*delte)
193 dent=entropy(indeold+1)-entropy(indeold)
194 cd write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
195 cd write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
197 sold=entropy(indeold)+(dent/delte)*deix
199 write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
200 write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,
202 if (minim) call zapis(varia,etot)
204 C NACC is the counter for the accepted conformations of a given processor
206 C NACC_TOT counts the total number of accepted conformations
209 if (MyID.eq.MasterID) then
210 call receive_MCM_info
212 call send_MCM_info(2)
215 do iene=indminn,indmaxx
222 c----------------------------------------------------------------------------
228 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
229 & 'Beginning iteration #',it
230 C Initialize local counter.
231 ntrial=0 ! # of generated non-overlapping confs.
232 noverlap=0 ! # of overlapping confs.
234 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
236 C Retrieve the angles of previously accepted conformation
240 cd write (iout,'(a)') 'Old variables:'
241 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
242 call var_to_geom(nvar,varia)
248 C Decide whether to generate a random conformation or perturb the old one
249 RandOrPert=ran_number(0.0D0,1.0D0)
250 if (RandOrPert.gt.RanFract) then
252 & write (iout,'(a)') 'Perturbation-generated conformation.'
253 call perturb(error,lprint,MoveType,nbond,1.0D0)
255 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
256 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
257 & MoveType,' returned from PERTURB.'
264 nstart_grow=iran_num(3,nres)
266 & write (iout,'(2a,i3)') 'Random-generated conformation',
267 & ' - chain regrown from residue',nstart_grow
269 call gen_rand_conf(nstart_grow,nrestmp,*30)
271 call geom_to_var(nvar,varia)
272 cd write (iout,'(a)') 'New variables:'
273 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
275 if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)')
276 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
277 if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)')
278 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
279 call etotal(energia(0))
281 c call enerprint(energia(0))
282 c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
283 if (etot-elowest.gt.overlap_cut) then
284 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,
285 & ' Overlap detected in the current conf.; energy is',etot
289 if (noverlap.gt.maxoverlap) then
290 write (iout,'(a)') 'Too many overlapping confs.'
295 call minimize(etot,varia,iretcode,nfun)
296 cd write (iout,'(a)') 'Variables after minimization:'
297 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
298 call etotal(energia(0))
300 neneval=neneval+nfun+1
302 if (print_mc.gt.2) then
303 write (iout,'(a)') 'Total energies of trial conf:'
304 call enerprint(energia(0))
305 else if (print_mc.eq.1) then
306 write (iout,'(a,i6,a,1pe16.6)')
307 & 'Trial conformation:',ngen,' energy:',etot
309 C--------------------------------------------------------------------------
311 C--------------------------------------------------------------------------
314 & call accepting(etot,eold,scur,sold,varia,varold,
319 if (elowest.gt.etot) elowest=etot
320 if (ehighest.lt.etot) ehighest=etot
321 moves_acc(MoveType)=moves_acc(MoveType)+1
322 if (MoveType.eq.1) then
323 nbond_acc(nbond)=nbond_acc(nbond)+1
325 C Check against conformation repetitions.
326 irep=conf_comp(varia,etot)
327 #if defined(AIX) || defined(PGI) || defined(CRAY)
328 open (istat,file=statname,position='append')
330 open (istat,file=statname,access='append')
334 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
336 & przes,obr,non_conv)
338 call contact(.false.,ncont,icont,co)
339 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
341 & write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
342 & 'RMS deviation from the reference structure:',rms,
343 & ' % of native contacts:',frac*100,' contact order:',co
345 & write (istat,'(i5,11(1pe14.5))') it,
346 & (energia(print_order(i)),i=1,nprint_ene),etot,
348 elseif (print_stat) then
349 write (istat,'(i5,10(1pe14.5))') it,
350 & (energia(print_order(i)),i=1,nprint_ene),etot
354 & call statprint(nacc,nfun,iretcode,etot,elowest)
355 C Print internal coordinates.
356 if (print_int) call briefout(nacc,etot)
358 if (MyID.ne.MasterID) then
359 call recv_stop_sig(Kwita)
360 cd print *,'Processor:',MyID,' STOP=',Kwita
362 call send_MCM_info(2)
364 call send_MCM_info(1)
368 C Store the accepted conf. and its energy.
376 cd write (iout,*) 'Accepted conformation:'
377 cd write (iout,*) (rad2deg*varia(i),i=1,nphi)
378 if (minim) call zapis(varia,etot)
380 ener(i,nsave)=energia(i)
382 ener(n_ene+1,nsave)=etot
383 ener(n_ene+2,nsave)=frac
385 nminima(irep)=nminima(irep)+1.0D0
386 c print *,'irep=',irep,' nminima=',nminima(irep)
388 if (Kwita.eq.0) call recv_stop_sig(kwita)
393 if (MyID.eq.MasterID) then
394 call receive_MCM_info
395 if (nacc_tot.ge.maxacc) accepted=.true.
398 if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
399 C Take a conformation from the pool
404 write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
406 & 'Take conformation',ii,' from the pool energy=',epool(ii)
408 & write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
410 endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
414 if (MyID.eq.MasterID) then
415 call receive_MCM_info
417 if (Kwita.eq.0) call recv_stop_sig(kwita)
419 if (ovrtim()) WhatsUp=-1
420 cd write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
421 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
423 cd write (iout,*) 'not_done=',not_done
426 print *,'Processor',MyID,
427 & ' has received STOP signal =',Kwita,' in EntSamp.'
428 cd print *,'not_done=',not_done
429 if (Kwita.lt.-1) WhatsUp=Kwita
430 else if (nacc_tot.ge.maxacc) then
431 print *,'Processor',MyID,' calls send_stop_sig,',
432 & ' because a sufficient # of confs. have been collected.'
433 cd print *,'not_done=',not_done
434 call send_stop_sig(-1)
435 else if (WhatsUp.eq.-1) then
436 print *,'Processor',MyID,
437 & ' calls send_stop_sig because of timeout.'
438 cd print *,'not_done=',not_done
439 call send_stop_sig(-2)
444 C-----------------------------------------------------------------
445 C... Construct energy histogram & update entropy
446 C-----------------------------------------------------------------
450 write (iout,*) 'Processor',MyID,
451 & ' is broadcasting ERROR-STOP signal.'
452 write (*,*) 'Processor',MyID,
453 & ' is broadcasting ERROR-STOP signal.'
454 call send_stop_sig(-3)
458 if (MyID.eq.MasterID) then
459 c call receive_MCM_results
460 call receive_energies
463 if (esave(i).lt.elowest) elowest=esave(i)
464 if (esave(i).gt.ehighest) ehighest=esave(i)
466 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
467 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
468 & ' Highest energy',ehighest
469 if (isweep.eq.1 .and. .not.ent_read) then
472 write (iout,*) 'EMAX=',emax
474 indmaxx=(ehighest-emin)/delte
477 do i=-max_ene,max_ene
478 entropy(i)=(emin+i*delte)*betbol
482 indmin=(elowest-emin)/delte
483 indmax=(ehighest-emin)/delte
484 if (indmin.lt.indminn) indminn=indmin
485 if (indmax.gt.indmaxx) indmaxx=indmax
487 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
488 C Construct energy histogram
490 inde=(esave(i)-emin)/delte
491 nhist(inde)=nhist(inde)+nminima(i)
493 C Update entropy (density of states)
495 if (nhist(i).gt.0) then
496 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
500 Cd entropy(i)=1.0D+10
502 write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
503 & 'End of macroiteration',isweep
504 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
505 & ' Ehighest=',ehighest
506 write (iout,'(a)') 'Frequecies of minima'
508 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
510 write (iout,'(/a)') 'Energy histogram'
512 write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
514 write (iout,'(/a)') 'Entropy'
516 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
518 C-----------------------------------------------------------------
519 C... End of energy histogram construction
520 C-----------------------------------------------------------------
522 entropy(-max_ene-4)=dfloat(indminn)
523 entropy(-max_ene-3)=dfloat(indmaxx)
524 entropy(-max_ene-2)=emin
525 entropy(-max_ene-1)=emax
527 cd print *,entname,ientout
528 open (ientout,file=entname,status='unknown')
529 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
531 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
535 write (iout,'(a)') 'Frequecies of minima'
537 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
539 c call send_MCM_results
541 call receive_MCM_update
542 indminn=entropy(-max_ene-4)
543 indmaxx=entropy(-max_ene-3)
544 emin=entropy(-max_ene-2)
545 emax=entropy(-max_ene-1)
546 write (iout,*) 'Received from master:'
547 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
548 & ' emin=',emin,' emax=',emax
549 write (iout,'(/a)') 'Entropy'
551 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
554 if (WhatsUp.lt.-1) return
556 if (ovrtim() .or. WhatsUp.lt.0) return
559 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
560 call statprint(nacc,nfun,iretcode,etot,elowest)
562 & 'Statistics of multiple-bond motions. Total motions:'
563 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
564 write (iout,'(a)') 'Accepted motions:'
565 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
566 write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
567 write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
569 C---------------------------------------------------------------------------
571 C---------------------------------------------------------------------------
575 if (isweep.eq.nsweep .and. it.ge.maxacc)
576 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
579 c------------------------------------------------------------------------------
580 subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
581 implicit real*8 (a-h,o-z)
585 include 'COMMON.IOUNITS'
588 include 'COMMON.INFO'
591 double precision ecur,eold,xx,ran_number,bol
592 double precision x(maxvar),xold(maxvar)
593 double precision tole /1.0D-1/, tola /5.0D0/
595 C Check if the conformation is similar.
596 cd write (iout,*) 'Enter ACCEPTING'
597 cd write (iout,*) 'Old PHI angles:'
598 cd write (iout,*) (rad2deg*xold(i),i=1,nphi)
599 cd write (iout,*) 'Current angles'
600 cd write (iout,*) (rad2deg*x(i),i=1,nphi)
601 cd ddif=dif_ang(nphi,x,xold)
602 cd write (iout,*) 'Angle norm:',ddif
603 cd write (iout,*) 'ecur=',ecur,' emax=',emax
604 if (ecur.gt.emax) then
607 & write (iout,'(a)') 'Conformation rejected as too high in energy'
609 else if (dabs(ecur-eold).lt.tole .and.
610 & dif_ang(nphi,x,xold).lt.tola) then
613 & write (iout,'(a)') 'Conformation rejected as too similar'
616 C Else evaluate the entropy of the conf and compare it with that of the previous
618 indecur=(ecur-emin)/delte
619 if (iabs(indecur).gt.max_ene) then
620 write (iout,'(a,2i5)')
621 & 'Accepting: Index out of range:',indecur
623 else if (indecur.eq.indmaxx) then
624 scur=entropy(indecur)
625 if (print_mc.gt.0) write (iout,*)'Energy boundary reached',
626 & indmaxx,indecur,entropy(indecur)
628 deix=ecur-(emin+indecur*delte)
629 dent=entropy(indecur+1)-entropy(indecur)
630 scur=entropy(indecur)+(dent/delte)*deix
632 cd print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
633 cd & ' scur=',scur,' eold=',eold,' sold=',sold
634 cd print *,'deix=',deix,' dent=',dent,' delte=',delte
635 if (print_mc.gt.1) then
636 write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
637 write(iout,*)'eold=',eold,' sold=',sold
639 if (scur.le.sold) then
642 C Else carry out acceptance test
643 xx=ran_number(0.0D0,1.0D0)
645 if (xxh.gt.50.0D0) then
652 if (print_mc.gt.0) write (iout,'(a)')
653 & 'Conformation accepted.'
656 if (print_mc.gt.0) write (iout,'(a)')
657 & 'Conformation rejected.'
662 c-----------------------------------------------------------------------------
664 implicit real*8 (a-h,o-z)
666 include 'COMMON.IOUNITS'
671 double precision varia(maxvar)
672 print '(a)','Call READ_POOL'
675 read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
676 if (epool(npool).eq.0.0D0) goto 10
677 call read_angles(intin,*10)
678 call geom_to_var(nvar,xpool(1,npool))
682 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
683 if (print_mc.gt.2) then
685 write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',
687 write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
689 endif ! (print_mc.gt.2)