2 C Does Boltzmann and entropic sampling without energy minimization
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,nbins
22 double precision RandOrPert
23 double precision varia(maxvar),elowest,elowest1,
24 & ehighest,ehighest1,eold
25 double precision przes(3),obr(3,3)
26 double precision varold(maxvar)
28 integer moves1(-1:MaxMoveType+1,0:MaxProcs-1),
29 & moves_acc1(-1:MaxMoveType+1,0:MaxProcs-1)
31 double precision etot_temp,etot_all(0:MaxProcs)
32 external d_vadd,d_vmin,d_vmax
33 double precision entropy1(-max_ene:max_ene),
34 & nhist1(-max_ene:max_ene)
35 integer nbond_move1(maxres*(MaxProcs+1)),
36 & nbond_acc1(maxres*(MaxProcs+1)),itemp(2)
38 double precision var_lowest(maxvar)
39 double precision energia(0:n_ene),energia_ave(0:n_ene)
41 write(iout,'(a,i8,2x,a,f10.5)')
42 & 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
43 open (istat,file=statname)
47 facee=1.0D0/(maxacc*delte)
48 C Number of bins in energy histogram
50 write (iout,*) 'NBINS=',nbins
52 C Read entropy from previous simulations.
54 read (ientin,*) indminn,indmaxx,emin,emax
55 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
60 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
63 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
64 & ' emin=',emin,' emax=',emax
65 write (iout,'(/a)') 'Initial entropy'
67 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
70 C Read the pool of conformations
74 C----------------------------------------------------------------------------
75 C Entropy-sampling simulations with continually updated entropy;
76 C set NSWEEP=1 for Boltzmann sampling.
77 C Loop thru simulations
78 C----------------------------------------------------------------------------
81 C Initialize the IFINISH array.
88 c---------------------------------------------------------------------------
89 C Initialize counters.
90 c---------------------------------------------------------------------------
91 C Total number of generated confs.
93 C Total number of moves. In general this won't be equal to the number of
94 C attempted moves, because we may want to reject some "bad" confs just by
97 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
103 C Initialize total and accepted number of moves of various kind.
108 C Total number of energy evaluations.
111 C----------------------------------------------------------------------------
112 C Take a conformation from the pool
113 C----------------------------------------------------------------------------
115 write (iout,*) 'emin=',emin,' emax=',emax
121 write (iout,*) 'Took conformation',ii,' from the pool energy=',
123 call var_to_geom(nvar,varia)
124 C Print internal coordinates of the initial conformation
126 else if (isweep.gt.1) then
127 if (eold.lt.emax) then
133 varia(i)=var_lowest(i)
136 call var_to_geom(nvar,varia)
138 C----------------------------------------------------------------------------
139 C Compute and print initial energies.
140 C----------------------------------------------------------------------------
144 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
145 write (iout,'(/80(1h*)/a)') 'Initial energies:'
147 call geom_to_var(nvar,varia)
148 call etotal(energia(0))
150 call enerprint(energia(0))
152 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,
155 call contact(.false.,ncont,icont,co)
156 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
157 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
158 & 'RMS deviation from the reference structure:',rms,
159 & ' % of native contacts:',frac*100,' contact order',co
160 write (istat,'(i10,16(1pe14.5))') 0,
161 & (energia(print_order(i)),i=1,nprint_ene),
164 write (istat,'(i10,14(1pe14.5))') 0,
165 & (energia(print_order(i)),i=1,nprint_ene),etot
169 if (.not. ent_read) then
170 C Initialize the entropy array
172 C Collect total energies from other processors.
175 call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
176 if (MyID.eq.MasterID) then
177 C Get the lowest and the highest energy.
178 print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),
179 & ' emin=',emin,' emax=',emax
183 if (emin.gt.etot_all(i)) emin=etot_all(i)
184 if (emax.lt.etot_all(i)) emax=etot_all(i)
187 endif ! MyID.eq.MasterID
190 print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
191 call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
192 print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
193 if (MyID.ne.MasterID) then
194 print *,'Processor:',MyID,etot_all(1),etot_all(2),
195 & etot_all(1),etot_all(2)
198 endif ! MyID.ne.MasterID
199 write (iout,*) 'After MP_GATHER etot_temp=',
200 & etot_temp,' emin=',emin
208 C Multicanonical sampling - start from Boltzmann distribution
209 do i=-max_ene,max_ene
210 entropy(i)=(emin+i*delte)*betbol
213 C Entropic sampling - start from uniform distribution of the density of states
214 do i=-max_ene,max_ene
218 write (iout,'(/a)') 'Initial entropy'
220 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
222 if (isweep.eq.1) then
226 indmaxx=indminn+nbins
229 endif ! .not. ent_read
231 call recv_stop_sig(Kwita)
232 if (whatsup.eq.1) then
233 call send_stop_sig(-2)
235 else if (whatsup.le.-2) then
237 else if (whatsup.eq.2) then
245 write (iout,'(/80(1h*)/20x,a/80(1h*))')
246 & 'Enter Monte Carlo procedure.'
248 call briefout(0,etot)
253 call entropia(eold,sold,indeold)
254 C NACC is the counter for the accepted conformations of a given processor
256 C NACC_TOT counts the total number of accepted conformations
259 c----------------------------------------------------------------------------
260 C Zero out average energies
264 C Initialize energy histogram
265 do i=-max_ene,max_ene
268 C Zero out iteration counter.
273 C Begin MC iteration loop.
276 C Initialize local counter.
277 ntrial=0 ! # of generated non-overlapping confs.
278 noverlap=0 ! # of overlapping confs.
280 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
282 C Retrieve the angles of previously accepted conformation
286 call var_to_geom(nvar,varia)
292 C Decide whether to take a conformation from the pool or generate/perturb one
294 from_pool=ran_number(0.0D0,1.0D0)
295 if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
296 C Throw a dice to choose the conformation from the pool
301 call var_to_geom(nvar,varia)
304 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
305 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
306 & write (iout,'(a,i3,a,f10.5)')
307 & 'Try conformation',ii,' from the pool energy=',epool(ii)
309 moves(-1)=moves(-1)+1
311 C Decide whether to generate a random conformation or perturb the old one
312 RandOrPert=ran_number(0.0D0,1.0D0)
313 if (RandOrPert.gt.RanFract) then
314 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
315 & write (iout,'(a)') 'Perturbation-generated conformation.'
316 call perturb(error,lprint,MoveType,nbond,0.1D0)
318 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
319 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
320 & MoveType,' returned from PERTURB.'
327 nstart_grow=iran_num(3,nres)
328 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
329 & write (iout,'(2a,i3)') 'Random-generated conformation',
330 & ' - chain regrown from residue',nstart_grow
332 call gen_rand_conf(nstart_grow,nrestmp,*30)
334 call geom_to_var(nvar,varia)
336 Cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
338 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
339 & write (iout,'(a,i5,a,i10,a,i10)')
340 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
341 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
342 & write (*,'(a,i5,a,i10,a,i10)')
343 & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
344 call etotal(energia(0))
347 cd call enerprint(energia(0))
348 cd write(iout,*)'it=',it,' etot=',etot
349 if (etot-elowest.gt.overlap_cut) then
350 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
351 & write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,
352 & ' Overlap detected in the current conf.; energy is',etot
355 if (noverlap.gt.maxoverlap) then
356 write (iout,'(a)') 'Too many overlapping confs.'
360 C--------------------------------------------------------------------------
362 C--------------------------------------------------------------------------
365 & call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
369 if (elowest.gt.etot) then
372 var_lowest(i)=varia(i)
375 if (ehighest.lt.etot) ehighest=etot
376 moves_acc(MoveType)=moves_acc(MoveType)+1
377 if (MoveType.eq.1) then
378 nbond_acc(nbond)=nbond_acc(nbond)+1
380 C Compare with reference structure.
382 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),
383 & nsup,przes,obr,non_conv)
385 call contact(.false.,ncont,icont,co)
386 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
389 C Periodically save average energies and confs.
392 energia_ave(i)=energia_ave(i)+energia(i)
394 moves(MaxMoveType+1)=nmove
395 moves_acc(MaxMoveType+1)=nacc
396 IF ((it/save_frequency)*save_frequency.eq.it) THEN
398 energia_ave(i)=energia_ave(i)/save_frequency
400 etot_ave=energia_ave(0)
402 C open (istat,file=statname,position='append')
404 C open (istat,file=statname,access='append')
407 & write (iout,'(80(1h*)/20x,a,i20)')
409 if (refstr .and. print_mc.gt.0) then
410 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
411 & 'RMS deviation from the reference structure:',rms,
412 & ' % of native contacts:',frac*100,' contact order:',co
416 write (istat,'(i10,10(1pe14.5))') it,
417 & (energia_ave(print_order(i)),i=1,nprint_ene),
418 & etot_ave,rms_ave,frac_ave
420 write (istat,'(i10,10(1pe14.5))') it,
421 & (energia_ave(print_order(i)),i=1,nprint_ene),
427 & call statprint(nacc,nfun,iretcode,etot,elowest)
428 C Print internal coordinates.
429 if (print_int) call briefout(nacc,etot)
433 ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
435 inde=icialosc((etot-emin)/delte)
436 nhist(inde)=nhist(inde)+1.0D0
438 if ( (it/message_frequency)*message_frequency.eq.it
439 & .and. (MyID.ne.MasterID) ) then
440 call recv_stop_sig(Kwita)
441 call send_MCM_info(message_frequency)
444 C Store the accepted conf. and its energy.
451 if (Kwita.eq.0) call recv_stop_sig(kwita)
456 if (MyID.eq.MasterID .and.
457 & (it/message_frequency)*message_frequency.eq.it) then
459 if (nacc_tot.ge.maxacc) accepted=.true.
462 C if ((ntrial.gt.maxtrial_iter
463 C & .or. (it/pool_read_freq)*pool_read_freq.eq.it)
464 C & .and. npool.gt.0) then
465 C Take a conformation from the pool
466 C ii=iran_num(1,npool)
468 C varold(i)=xpool(i,ii)
470 C if (ntrial.gt.maxtrial_iter)
471 C & write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
473 C & 'Take conformation',ii,' from the pool energy=',epool(ii)
475 C & write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
478 C call entropia(eold,sold,indeold)
480 C endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
484 if (MyID.eq.MasterID .and.
485 & (it/message_frequency)*message_frequency.eq.it) then
488 if (Kwita.eq.0) call recv_stop_sig(kwita)
490 if (ovrtim()) WhatsUp=-1
491 cd write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
492 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
494 cd write (iout,*) 'not_done=',not_done
497 print *,'Processor',MyID,
498 & ' has received STOP signal =',Kwita,' in EntSamp.'
499 cd print *,'not_done=',not_done
500 if (Kwita.lt.-1) WhatsUp=Kwita
501 if (MyID.ne.MasterID) call send_MCM_info(-1)
502 else if (nacc_tot.ge.maxacc) then
503 print *,'Processor',MyID,' calls send_stop_sig,',
504 & ' because a sufficient # of confs. have been collected.'
505 cd print *,'not_done=',not_done
506 call send_stop_sig(-1)
507 if (MyID.ne.MasterID) call send_MCM_info(-1)
508 else if (WhatsUp.eq.-1) then
509 print *,'Processor',MyID,
510 & ' calls send_stop_sig because of timeout.'
511 cd print *,'not_done=',not_done
512 call send_stop_sig(-2)
513 if (MyID.ne.MasterID) call send_MCM_info(-1)
518 C-----------------------------------------------------------------
519 C... Construct energy histogram & update entropy
520 C-----------------------------------------------------------------
524 write (iout,*) 'Processor',MyID,
525 & ' is broadcasting ERROR-STOP signal.'
526 write (*,*) 'Processor',MyID,
527 & ' is broadcasting ERROR-STOP signal.'
528 call send_stop_sig(-3)
529 if (MyID.ne.MasterID) call send_MCM_info(-1)
532 write (iout,'(/a)') 'Energy histogram'
534 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
537 C Wait until every processor has sent complete MC info.
538 if (MyID.eq.MasterID) then
541 C write (*,*) 'The IFINISH array:'
542 C write (*,*) (ifinish(i),i=1,nctasks)
545 not_done=not_done.or.(ifinish(i).ge.0)
547 if (not_done) call receive_MC_info
550 C Make collective histogram from the work of all processors.
551 msglen=(2*max_ene+1)*8
553 & 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',
555 call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,
557 print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
558 do i=-max_ene,max_ene
561 C Collect min. and max. energy
563 &'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
564 call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
565 call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
566 print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
567 IF (MyID.eq.MasterID) THEN
571 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
572 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
573 & ' Highest energy',ehighest
574 indmin=icialosc((elowest-emin)/delte)
575 imdmax=icialosc((ehighest-emin)/delte)
576 if (indmin.lt.indminn) then
577 emax=emin+indmin*delte+e_up
581 if (.not.ent_read) ent_read=.true.
582 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
583 C Update entropy (density of states)
585 if (nhist(i).gt.0) then
586 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
589 write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
590 & 'End of macroiteration',isweep
591 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
592 & ' Ehighest=',ehighest
593 write (iout,'(/a)') 'Energy histogram'
595 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
597 write (iout,'(/a)') 'Entropy'
599 write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
601 C-----------------------------------------------------------------
602 C... End of energy histogram construction
603 C-----------------------------------------------------------------
606 if (.not. ent_read) ent_read=.true.
607 ENDIF ! MyID .eq. MaterID
608 if (MyID.eq.MasterID) then
612 print *,'before mp_bcast processor',MyID,' indminn=',indminn,
613 & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
614 call mp_bcast(itemp(1),8,MasterID,cgGroupID)
615 call mp_bcast(emax,8,MasterID,cgGroupID)
616 print *,'after mp_bcast processor',MyID,' indminn=',indminn,
617 & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
618 if (MyID .ne. MasterID) then
622 msglen=(indmaxx-indminn+1)*8
623 print *,'processor',MyID,' calling mp_bcast msglen=',msglen,
624 & ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
625 call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
626 IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
627 open (ientout,file=entname,status='unknown')
628 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
630 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
634 write (iout,*) 'Received from master:'
635 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
636 & ' emin=',emin,' emax=',emax
637 write (iout,'(/a)') 'Entropy'
639 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
641 ENDIF ! MyID.eq.MasterID
642 print *,'Processor',MyID,' calls MP_GATHER'
643 call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,
645 call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,
647 print *,'Processor',MyID,' MP_GATHER call accomplished'
648 if (MyID.eq.MasterID) then
650 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
651 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
653 & 'Statistics of multiple-bond motions. Total motions:'
654 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
655 write (iout,'(a)') 'Accepted motions:'
656 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
659 & 'Statistics of multi-bond moves of respective processors:'
663 nbond_move(i)=nbond_move(i)+nbond_move1(ind)
664 nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
668 write (iout,*) 'Processor',iproc,' nbond_move:',
669 & (nbond_move1(iproc*nbm+i),i=1,Nbm),
670 & ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
673 call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,
675 call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),
676 & MasterID,cgGroupID)
677 if (MyID.eq.MasterID) then
679 do i=-1,MaxMoveType+1
680 moves(i)=moves(i)+moves1(i,iproc)
681 moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
689 write (iout,*) 'Processor',iproc,' moves',
690 & (moves1(i,iproc),i=0,MaxMoveType+1),
691 & ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
695 open (ientout,file=entname,status='unknown')
696 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
698 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
702 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
703 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
705 & 'Statistics of multiple-bond motions. Total motions:'
706 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
707 write (iout,'(a)') 'Accepted motions:'
708 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
709 if (ovrtim() .or. WhatsUp.lt.0) return
711 C---------------------------------------------------------------------------
713 C---------------------------------------------------------------------------
717 if (isweep.eq.nsweep .and. it.ge.maxacc)
718 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
721 c------------------------------------------------------------------------------
722 subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
723 implicit real*8 (a-h,o-z)
727 include 'COMMON.IOUNITS'
730 include 'COMMON.INFO'
733 double precision ecur,eold,xx,ran_number,bol
734 double precision x(maxvar),xold(maxvar)
736 C Check if the conformation is similar.
737 cd write (iout,*) 'Enter ACCEPTING'
738 cd write (iout,*) 'Old PHI angles:'
739 cd write (iout,*) (rad2deg*xold(i),i=1,nphi)
740 cd write (iout,*) 'Current angles'
741 cd write (iout,*) (rad2deg*x(i),i=1,nphi)
742 cd ddif=dif_ang(nphi,x,xold)
743 cd write (iout,*) 'Angle norm:',ddif
744 cd write (iout,*) 'ecur=',ecur,' emax=',emax
745 if (ecur.gt.emax) then
747 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
748 & write (iout,'(a)') 'Conformation rejected as too high in energy'
751 C Else evaluate the entropy of the conf and compare it with that of the previous
753 call entropia(ecur,scur,indecur)
754 cd print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
755 cd & ' scur=',scur,' eold=',eold,' sold=',sold
756 cd print *,'deix=',deix,' dent=',dent,' delte=',delte
757 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
758 write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,
760 write(iout,*)'eold=',eold,' sold=',sold
762 if (scur.le.sold) then
765 C Else carry out acceptance test
766 xx=ran_number(0.0D0,1.0D0)
768 if (xxh.gt.50.0D0) then
775 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
776 & write (iout,'(a)') 'Conformation accepted.'
779 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
780 & write (iout,'(a)') 'Conformation rejected.'
785 c--------------------------------------------------------------------------
786 integer function icialosc(x)
795 c--------------------------------------------------------------------------
796 subroutine entropia(ecur,scur,indecur)
797 implicit real*8 (a-h,o-z)
801 include 'COMMON.IOUNITS'
802 double precision ecur,scur
804 indecur=icialosc((ecur-emin)/delte)
805 if (iabs(indecur).gt.max_ene) then
806 if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)')
807 & 'Accepting: Index out of range:',indecur
809 else if (indecur.ge.indmaxx) then
810 scur=entropy(indecur)
811 if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it)
812 & write (iout,*)'Energy boundary reached',
813 & indmaxx,indecur,entropy(indecur)
815 deix=ecur-(emin+indecur*delte)
816 dent=entropy(indecur+1)-entropy(indecur)
817 scur=entropy(indecur)+(dent/delte)*deix