2 implicit real*8 (a-h,o-z)
4 include 'COMMON.IOUNITS'
6 include 'COMMON.CONTROL'
7 include 'COMMON.INTERACT'
12 C Set up variables used in MC/MCM.
14 write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
15 write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,
16 & ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
17 write (iout,'(4(a,f8.1)/2(a,i3))')
18 & 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,
19 & ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
20 if (nwindow.gt.0) then
21 write (iout,'(a)') 'Perturbation windows:'
27 write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,
31 C Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
33 C Number of "end bonds".
36 print *,'koniecl=',koniecl
37 write (iout,'(a)') 'Probabilities of move types:'
38 write (*,'(a)') 'Probabilities of move types:'
40 write (iout,'(a,f10.5)') MovTypID(i),
41 & sumpro_type(i)-sumpro_type(i-1)
42 write (*,'(a,f10.5)') MovTypID(i),
43 & sumpro_type(i)-sumpro_type(i-1)
46 C Maximum length of N-bond segment to be moved
47 c nbm=nres-1-(2*koniecl-1)
48 if (nwindow.gt.0) then
51 if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
54 write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
61 sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
63 write (iout,'(a)') 'The SumPro_Bond array:'
64 write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
65 write (*,'(a)') 'The SumPro_Bond array:'
66 write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
67 C Maximum number of side chains moved simultaneously
68 c print *,'nnt=',nnt,' nct=',nct
71 if (itype(i).eq.10) ngly=ngly+1
75 MaxSideMove=min0((nct-nnt+1)/2,mmm)
77 c print *,'MaxSideMove=',MaxSideMove
78 C Max. number of generated confs (not used at present).
80 C Set initial temperature
82 betbol=1.0D0/(Rbol*Tcur)
83 write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,
85 write (iout,*) 'RanFract=',ranfract
88 c------------------------------------------------------------------------------
90 subroutine do_mcm(i_orig)
91 C Monte-Carlo-with-Minimization calculations - serial code.
92 implicit real*8 (a-h,o-z)
94 include 'COMMON.IOUNITS'
96 include 'COMMON.CHAIN'
98 include 'COMMON.CONTACTS'
99 include 'COMMON.CONTROL'
101 include 'COMMON.INTERACT'
102 include 'COMMON.CACHE'
103 crc include 'COMMON.DEFORM'
104 crc include 'COMMON.DEFORM1'
105 include 'COMMON.NAMES'
106 logical accepted,over,ovrtim,error,lprint,not_done,my_conf,
108 integer MoveType,nbond,conf_comp
109 integer ifeed(max_cache)
110 double precision varia(maxvar),varold(maxvar),elowest,eold,
112 double precision energia(0:n_ene)
113 double precision coord1(maxres,3)
116 C---------------------------------------------------------------------------
117 C Initialize counters.
118 C---------------------------------------------------------------------------
119 C Total number of generated confs.
121 C Total number of moves. In general this won't be equal to the number of
122 C attempted moves, because we may want to reject some "bad" confs just by
125 C Total number of temperature jumps.
127 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
133 C Initialize total and accepted number of moves of various kind.
138 C Total number of energy evaluations.
143 write (iout,*) 'RanFract=',RanFract
148 c----------------------------------------------------------------------------
149 C Compute and print initial energies.
150 c----------------------------------------------------------------------------
152 write (iout,'(/80(1h*)/a)') 'Initial energies:'
156 call etotal(energia(0))
158 C Minimize the energy of the first conformation.
160 call geom_to_var(nvar,varia)
161 ! write (iout,*) 'The VARIA array'
162 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
163 call minimize(etot,varia,iretcode,nfun)
164 call var_to_geom(nvar,varia)
166 write (iout,*) 'etot from MINIMIZE:',etot
167 ! write (iout,*) 'Tha VARIA array'
168 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
170 call etotal(energia(0))
172 call enerprint(energia(0))
175 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
178 call contact(.false.,ncont,icont,co)
179 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
180 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
181 & 'RMS deviation from the reference structure:',rms,
182 & ' % of native contacts:',frac*100,' contact order:',co
184 & write (istat,'(i5,17(1pe14.5))') 0,
185 & (energia(print_order(i)),i=1,nprint_ene),
188 if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,
189 & (energia(print_order(i)),i=1,nprint_ene),etot
191 if (print_stat) close(istat)
192 neneval=neneval+nfun+1
193 write (iout,'(/80(1h*)/20x,a/80(1h*))')
194 & 'Enter Monte Carlo procedure.'
197 call briefout(0,etot)
204 call zapis(varia,etot)
205 nacc=0 ! total # of accepted confs of the current processor.
206 nacc_tot=0 ! total # of accepted confs of all processors.
208 not_done = (iretcode.ne.11)
210 C----------------------------------------------------------------------------
212 c----------------------------------------------------------------------------
217 write (iout,'(80(1h*)/20x,a,i7)')
218 & 'Beginning iteration #',it
219 C Initialize local counter.
220 ntrial=0 ! # of generated non-overlapping confs.
222 do while (.not. accepted)
224 C Retrieve the angles of previously accepted conformation
225 noverlap=0 ! # of overlapping confs.
229 call var_to_geom(nvar,varia)
232 C Heat up the system, if necessary.
234 C If temperature cannot be further increased, stop.
239 cd write (iout,'(a)') 'Old variables:'
240 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
241 C Decide whether to generate a random conformation or perturb the old one
242 RandOrPert=ran_number(0.0D0,1.0D0)
243 if (RandOrPert.gt.RanFract) then
245 & write (iout,'(a)') 'Perturbation-generated conformation.'
246 call perturb(error,lprint,MoveType,nbond,1.0D0)
248 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
249 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
250 & MoveType,' returned from PERTURB.'
257 nstart_grow=iran_num(3,nres)
259 & write (iout,'(2a,i3)') 'Random-generated conformation',
260 & ' - chain regrown from residue',nstart_grow
261 call gen_rand_conf(nstart_grow,*30)
263 call geom_to_var(nvar,varia)
264 cd write (iout,'(a)') 'New variables:'
265 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
268 call etotal(energia(0))
270 c call enerprint(energia(0))
271 c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
272 if (etot-elowest.gt.overlap_cut) then
273 if(iprint.gt.1.or.etot.lt.1d20)
274 & write (iout,'(a,1pe14.5)')
275 & 'Overlap detected in the current conf.; energy is',etot
279 if (noverlap.gt.maxoverlap) then
280 write (iout,'(a)') 'Too many overlapping confs.'
285 call minimize(etot,varia,iretcode,nfun)
286 cd write (iout,*) 'etot from MINIMIZE:',etot
287 cd write (iout,'(a)') 'Variables after minimization:'
288 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
290 call etotal(energia(0))
292 neneval=neneval+nfun+2
294 c call enerprint(energia(0))
295 write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,
297 C--------------------------------------------------------------------------
298 C... Do Metropolis test
299 C--------------------------------------------------------------------------
303 if (WhatsUp.eq.0 .and. Kwita.eq.0) then
304 call metropolis(nvar,varia,varold,etot,eold,accepted,
305 & my_conf,EneLower,it)
307 write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
312 if (elowest.gt.etot) elowest=etot
313 moves_acc(MoveType)=moves_acc(MoveType)+1
314 if (MoveType.eq.1) then
315 nbond_acc(nbond)=nbond_acc(nbond)+1
317 C Check against conformation repetitions.
318 irepet=conf_comp(varia,etot)
320 #if defined(AIX) || defined(PGI) || defined(CRAY)
321 open (istat,file=statname,position='append')
323 open (istat,file=statname,access='append')
326 call statprint(nacc,nfun,iretcode,etot,elowest)
328 call var_to_geom(nvar,varia)
330 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
331 & nsup,przes,obr,non_conv)
333 call contact(.false.,ncont,icont,co)
334 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
335 write (iout,'(a,f8.3,a,f8.3)')
336 & 'RMS deviation from the reference structure:',rms,
337 & ' % of native contacts:',frac*100,' contact order',co
341 write (iout,*) 'Writing new conformation',nout
343 write (istat,'(i5,16(1pe14.5))') nout,
344 & (energia(print_order(i)),i=1,nprint_ene),
348 & write (istat,'(i5,17(1pe14.5))') nout,
349 & (energia(print_order(i)),i=1,nprint_ene),etot
351 if (print_stat) close(istat)
352 C Print internal coordinates.
353 if (print_int) call briefout(nout,etot)
354 C Accumulate the newly accepted conf in the coord1 array, if it is different
355 C from all confs that are already there.
356 call compare_s1(n_thr,max_thread2,etot,varia,ii,
357 & enetb1,coord1,rms_deform,.true.,iprint)
358 write (iout,*) 'After compare_ss: n_thr=',n_thr
359 if (ii.eq.1 .or. ii.eq.3) then
360 write (iout,'(8f10.4)')
361 & (rad2deg*coord1(i,n_thr),i=1,nvar)
364 write (iout,*) 'Conformation from cache, not written.'
367 if (nrepm.gt.maxrepm) then
368 write (iout,'(a)') 'Too many conformation repetitions.'
371 C Store the accepted conf. and its energy.
376 if (irepet.eq.0) call zapis(varia,etot)
377 C Lower the temperature, if necessary.
388 C Check for time limit.
389 if (ovrtim()) WhatsUp=-1
390 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
399 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
400 call statprint(nacc,nfun,iretcode,etot,elowest)
402 & 'Statistics of multiple-bond motions. Total motions:'
403 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
404 write (iout,'(a)') 'Accepted motions:'
405 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
407 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
413 c------------------------------------------------------------------------------
414 subroutine do_mcm(i_orig)
415 C Monte-Carlo-with-Minimization calculations - parallel code.
416 implicit real*8 (a-h,o-z)
419 include 'COMMON.IOUNITS'
421 include 'COMMON.CHAIN'
423 include 'COMMON.CONTACTS'
424 include 'COMMON.CONTROL'
426 include 'COMMON.INTERACT'
427 include 'COMMON.INFO'
428 include 'COMMON.CACHE'
429 crc include 'COMMON.DEFORM'
430 crc include 'COMMON.DEFORM1'
431 crc include 'COMMON.DEFORM2'
432 include 'COMMON.MINIM'
433 include 'COMMON.NAMES'
434 logical accepted,over,ovrtim,error,lprint,not_done,similar,
435 & enelower,non_conv,flag,finish
436 integer MoveType,nbond,conf_comp
437 double precision varia(maxvar),varold(maxvar),elowest,eold,
438 & x1(maxvar), varold1(maxvar), przes(3),obr(3,3)
439 integer iparentx(max_threadss2)
440 integer iparentx1(max_threadss2)
441 integer imtasks(150),imtasks_n
442 double precision energia(0:n_ene)
444 print *,'Master entered DO_MCM'
452 C---------------------------------------------------------------------------
453 C Initialize counters.
454 C---------------------------------------------------------------------------
455 C Total number of generated confs.
457 C Total number of moves. In general this won`t be equal to the number of
458 C attempted moves, because we may want to reject some "bad" confs just by
461 C Total number of temperature jumps.
463 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
469 C Initialize total and accepted number of moves of various kind.
474 C Total number of energy evaluations.
478 c write (iout,*) 'RanFract=',RanFract
481 c----------------------------------------------------------------------------
482 C Compute and print initial energies.
483 c----------------------------------------------------------------------------
485 write (iout,'(/80(1h*)/a)') 'Initial energies:'
488 call etotal(energia(0))
490 call enerprint(energia(0))
491 C Request energy computation from slave processors.
492 call geom_to_var(nvar,varia)
493 ! write (iout,*) 'The VARIA array'
494 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
495 call minimize(etot,varia,iretcode,nfun)
496 call var_to_geom(nvar,varia)
498 write (iout,*) 'etot from MINIMIZE:',etot
499 ! write (iout,*) 'Tha VARIA array'
500 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
503 if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))')
504 & 'Enter Monte Carlo procedure.'
505 if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
511 call zapis(varia,etot)
513 call var_to_geom(nvar,varia)
515 call etotal(energia(0))
516 if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
518 nacc=0 ! total # of accepted confs of the current processor.
519 nacc_tot=0 ! total # of accepted confs of all processors.
521 C----------------------------------------------------------------------------
523 c----------------------------------------------------------------------------
526 LOOP1:do while (not_done)
528 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
529 & 'Beginning iteration #',it
530 C Initialize local counter.
531 ntrial=0 ! # of generated non-overlapping confs.
532 noverlap=0 ! # of overlapping confs.
534 LOOP2:do while (.not. accepted)
536 LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
538 if(imtasks(i).eq.0) then
543 C Retrieve the angles of previously accepted conformation
547 call var_to_geom(nvar,varia)
550 C Heat up the system, if necessary.
552 C If temperature cannot be further increased, stop.
558 c write (iout,'(a)') 'Old variables:'
559 c write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
560 C Decide whether to generate a random conformation or perturb the old one
561 RandOrPert=ran_number(0.0D0,1.0D0)
562 if (RandOrPert.gt.RanFract) then
564 & write (iout,'(a)') 'Perturbation-generated conformation.'
565 call perturb(error,lprint,MoveType,nbond,1.0D0)
566 c print *,'after perturb',error,finish
567 if (error) finish = .true.
568 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
569 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
570 & MoveType,' returned from PERTURB.'
572 write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',
573 & MoveType,' returned from PERTURB.'
579 nstart_grow=iran_num(3,nres)
581 & write (iout,'(2a,i3)') 'Random-generated conformation',
582 & ' - chain regrown from residue',nstart_grow
583 call gen_rand_conf(nstart_grow,*30)
585 call geom_to_var(nvar,varia)
587 c print *,'finish=',finish
588 if (etot-elowest.gt.overlap_cut) then
589 if (print_mc.gt.1) write (iout,'(a,1pe14.5)')
590 & 'Overlap detected in the current conf.; energy is',etot
591 if(iprint.gt.1.or.etot.lt.1d19) print *,
592 & 'Overlap detected in the current conf.; energy is',etot
596 if (noverlap.gt.maxoverlap) then
597 write (iout,*) 'Too many overlapping confs.',
598 & ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
601 else if (.not. finish) then
602 C Distribute tasks to processors
603 c print *,'Master sending order'
604 call MPI_SEND(12, 1, MPI_INTEGER, is, tag,
606 c write (iout,*) '12: tag=',tag
607 c print *,'Master sent order to processor',is
608 call MPI_SEND(it, 1, MPI_INTEGER, is, tag,
610 c write (iout,*) 'it: tag=',tag
611 call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,
613 c write (iout,*) 'eold: tag=',tag
614 call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION,
617 c write (iout,*) 'varia: tag=',tag
618 call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION,
621 c write (iout,*) 'varold: tag=',tag
628 imtasks_n=imtasks_n+1
634 LOOP_RECV:do while(.not.flag)
636 call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
638 call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,
639 & CG_COMM, status, ierr)
640 call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,
641 & CG_COMM, status, ierr)
642 call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,
643 & CG_COMM, status, ierr)
644 call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,
645 & CG_COMM, status, ierr)
646 call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is,
647 & tag, CG_COMM, status, ierr)
648 call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,
649 & CG_COMM, status, ierr)
650 call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,
651 & CG_COMM, status, ierr)
652 call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,
653 & CG_COMM, status, ierr)
654 i_grnum_d=i_grnum_d+ii_grnum_d
655 i_ennum_d=i_ennum_d+ii_ennum_d
656 neneval = neneval+ii_ennum_d
657 i_hesnum_d=i_hesnum_d+ii_hesnum_d
658 i_minimiz=i_minimiz+1
660 imtasks_n=imtasks_n-1
666 if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)')
667 & 'From Worker #',is,' iitt',iitt,
668 & ' Conformation:',ngen,' energy:',etot
669 C--------------------------------------------------------------------------
670 C... Do Metropolis test
671 C--------------------------------------------------------------------------
672 call metropolis(nvar,varia,varold1,etot,eold1,accepted,
674 if(iitt.ne.it.and..not.similar) then
675 call metropolis(nvar,varia,varold,etot,eold,accepted,
679 if(etot.lt.eneglobal)eneglobal=etot
680 c if(mod(it,100).eq.0)
681 write(iout,*)'CHUJOJEB ',neneval,eneglobal
683 C Write the accepted conformation.
686 call var_to_geom(nvar,varia)
688 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
689 & nsup,przes,obr,non_conv)
691 call contact(.false.,ncont,icont,co)
692 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
693 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
694 & 'RMS deviation from the reference structure:',rms,
695 & ' % of native contacts:',frac*100,' contact order:',co
698 & write (iout,*) 'Writing new conformation',nout
700 call var_to_geom(nvar,varia)
701 #if defined(AIX) || defined(PGI) || defined(CRAY)
702 open (istat,file=statname,position='append')
704 open (istat,file=statname,access='append')
707 write (istat,'(i5,16(1pe14.5))') nout,
708 & (energia(print_order(i)),i=1,nprint_ene),
711 write (istat,'(i5,16(1pe14.5))') nout,
712 & (energia(print_order(i)),i=1,nprint_ene),etot
716 C Print internal coordinates.
717 if (print_int) call briefout(nout,etot)
720 if (elowest.gt.etot) elowest=etot
721 moves_acc(MoveType)=moves_acc(MoveType)+1
722 if (MoveType.eq.1) then
723 nbond_acc(nbond)=nbond_acc(nbond)+1
725 C Check against conformation repetitions.
726 irepet=conf_comp(varia,etot)
727 if (nrepm.gt.maxrepm) then
729 & write (iout,'(a)') 'Too many conformation repetitions.'
732 C Store the accepted conf. and its energy.
737 if (irepet.eq.0) call zapis(varia,etot)
738 C Lower the temperature, if necessary.
744 if(finish.and.imtasks_n.eq.0)exit LOOP2
745 enddo LOOP2 ! accepted
746 C Check for time limit.
747 not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
748 if(.not.not_done .or. finish) then
749 if(imtasks_n.gt.0) then
756 enddo LOOP1 ! not_done
758 if (print_mc.gt.0) then
759 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
760 call statprint(nacc,nfun,iretcode,etot,elowest)
762 & 'Statistics of multiple-bond motions. Total motions:'
763 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
764 write (iout,'(a)') 'Accepted motions:'
765 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
767 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
775 call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
780 c------------------------------------------------------------------------------
781 subroutine execute_slave(nodeinfo,iprint)
782 implicit real*8 (a-h,o-z)
785 include 'COMMON.TIME1'
786 include 'COMMON.IOUNITS'
787 crc include 'COMMON.DEFORM'
788 crc include 'COMMON.DEFORM1'
789 crc include 'COMMON.DEFORM2'
790 include 'COMMON.LOCAL'
792 include 'COMMON.INFO'
793 include 'COMMON.MINIM'
794 character*10 nodeinfo
795 double precision x(maxvar),x1(maxvar)
797 c print *,'Processor:',MyID,' Entering execute_slave'
799 c call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
802 1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
803 & CG_COMM, status, ierr)
804 c write(iout,*)'12: tag=',tag
805 if(iprint.ge.2)print *, MyID,' recv order ',i_switch
806 if (i_switch.eq.12) then
810 call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
811 & CG_COMM, status, ierr)
812 c write(iout,*)'12: tag=',tag
813 call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
814 & CG_COMM, status, ierr)
815 c write(iout,*)'ener: tag=',tag
816 call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
817 & CG_COMM, status, ierr)
818 c write(iout,*)'x: tag=',tag
819 call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
820 & CG_COMM, status, ierr)
821 c write(iout,*)'x1: tag=',tag
827 c print *,'calling minimize'
828 call minimize(energyx,x,iretcode,nfun)
830 & write(iout,100)'minimized energy = ',energyx,
831 & ' # funeval:',nfun,' iret ',iretcode
832 write(*,100)'minimized energy = ',energyx,
833 & ' # funeval:',nfun,' iret ',iretcode
834 100 format(a20,f10.5,a12,i5,a6,i2)
835 if(iretcode.eq.10) then
838 & write(iout,*)' ... not converged - trying again ',iminrep
839 call minimize(energyx,x,iretcode,nfun)
841 & write(iout,*)'minimized energy = ',energyx,
842 & ' # funeval:',nfun,' iret ',iretcode
843 if(iretcode.ne.10)go to 812
845 if(iretcode.eq.10) then
847 & write(iout,*)' ... not converged again - giving up'
852 c print *,'Sending results'
853 call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
855 call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
857 call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
859 call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
861 call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
863 call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
865 call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
867 call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
869 c print *,'End sending'
876 c------------------------------------------------------------------------------
877 subroutine statprint(it,nfun,iretcode,etot,elowest)
878 implicit real*8 (a-h,o-z)
880 include 'COMMON.IOUNITS'
881 include 'COMMON.CONTROL'
885 & '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
886 & 'Finished iteration #',it,' energy is',etot,
887 & ' lowest energy:',elowest,
888 & 'SUMSL return code:',iretcode,
889 & ' # of energy evaluations:',neneval,
890 & '# of temperature jumps:',ntherm,
891 & ' # of minima repetitions:',nrepm
893 write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
894 & 'Finished iteration #',it,' energy is',etot,
895 & ' lowest energy:',elowest
898 & 'Kind of move ',' total',' accepted',
900 write (iout,'(58(1h-))')
902 if (moves(i).eq.0) then
905 fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
907 write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
910 write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,
911 & dfloat(nacc_tot)/dfloat(nmove)
912 write (iout,'(58(1h-))')
913 write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
916 c------------------------------------------------------------------------------
917 subroutine heat(over)
918 implicit real*8 (a-h,o-z)
921 include 'COMMON.IOUNITS'
923 C Check if there`s a need to increase temperature.
924 if (ntrial.gt.maxtrial) then
925 if (NstepH.gt.0) then
926 if (dabs(Tcur-TMax).lt.1.0D-7) then
928 & write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))')
929 & 'Upper limit of temperature reached. Terminating.'
934 if (Tcur.gt.Tmax) Tcur=Tmax
935 betbol=1.0D0/(Rbol*Tcur)
937 & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
938 & 'System heated up to ',Tcur,' K; BetBol:',betbol
946 & 'Maximum number of trials in a single MCM iteration exceeded.'
955 c------------------------------------------------------------------------------
957 implicit real*8 (a-h,o-z)
960 include 'COMMON.IOUNITS'
961 if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
963 if (Tcur.lt.Tmin) Tcur=Tmin
964 betbol=1.0D0/(Rbol*Tcur)
966 & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
967 & 'System cooled down up to ',Tcur,' K; BetBol:',betbol
971 C---------------------------------------------------------------------------
972 subroutine zapis(varia,etot)
973 implicit real*8 (a-h,o-z)
977 include 'COMMON.INFO'
982 include 'COMMON.IOUNITS'
983 integer itemp(maxsave)
984 double precision varia(maxvar)
988 write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
989 & ' MaxSave=',MaxSave
990 write (iout,'(a)') 'Current energy and conformation:'
991 write (iout,'(1pe14.5)') etot
992 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
994 C Shift the contents of the esave and varsave arrays if filled up.
995 call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
996 & etot,varia,esave,varsave)
998 write (iout,'(a)') 'Energies and the VarSave array.'
1000 write (iout,'(i5,1pe14.5)') i,esave(i)
1001 write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
1006 C---------------------------------------------------------------------------
1007 subroutine perturb(error,lprint,MoveType,nbond,max_phi)
1008 implicit real*8 (a-h,o-z)
1009 include 'DIMENSIONS'
1010 parameter (MMaxSideMove=100)
1011 include 'COMMON.MCM'
1012 include 'COMMON.CHAIN'
1013 include 'COMMON.INTERACT'
1014 include 'COMMON.VAR'
1015 include 'COMMON.GEO'
1016 include 'COMMON.NAMES'
1017 include 'COMMON.IOUNITS'
1018 crc include 'COMMON.DEFORM1'
1019 logical error,lprint,fail
1020 integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
1021 double precision max_phi
1022 double precision psi,gen_psi
1029 C Perturb the conformation according to a randomly selected move.
1030 call SelectMove(MoveType)
1031 c write (iout,*) 'MoveType=',MoveType
1033 goto (100,200,300,400,500) MoveType
1034 C------------------------------------------------------------------------------
1035 C Backbone N-bond move.
1036 C Select the number of bonds (length of the segment to perturb).
1038 if (itrial.gt.1000) then
1039 write (iout,'(a)') 'Too many attempts at multiple-bond move.'
1043 bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
1044 c print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
1045 c & ' Bond_prob=',Bond_Prob
1047 c print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
1048 if (bond_prob.ge.sumpro_bond(i) .and.
1049 & bond_prob.le.sumpro_bond(i+1)) then
1054 write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
1055 & ' to move out of range.'
1059 if (nwindow.gt.0) then
1060 C Select the first residue to perturb
1061 iwindow=iran_num(1,nwindow)
1062 print *,'iwindow=',iwindow
1064 do while (winlen(iwindow).lt.nbond)
1065 iwindow=iran_num(1,nwindow)
1067 if (iiwin.gt.1000) then
1068 write (iout,'(a)') 'Cannot select moveable residues.'
1073 nstart=iran_num(winstart(iwindow),winend(iwindow))
1075 nstart = iran_num(koniecl+2,nres-nbond-koniecl)
1076 cd print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
1077 cd & ' nstart=',nstart
1080 if (psi.eq.0.0) then
1084 if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
1085 & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
1086 cd print *,'nstart=',nstart
1087 call bond_move(nbond,nstart,psi,.false.,error)
1090 & 'Could not define reference system in bond_move, ',
1091 & 'choosing ahother segment.'
1095 nbond_move(nbond)=nbond_move(nbond)+1
1099 C------------------------------------------------------------------------------
1100 C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
1104 c end_select=iran_num(1,2*koniecl)
1105 c if (end_select.gt.koniecl) then
1106 c end_select=nphi-(end_select-koniecl)
1108 c end_select=koniecl+3
1110 c if (nwindow.gt.0) then
1111 c iwin=iran_num(1,nwindow)
1112 c i1=max0(4,winstart(iwin))
1113 c i2=min0(winend(imin)+2,nres)
1114 c end_select=iran_num(i1,i2)
1116 c iselect = iran_num(1,nmov_var)
1119 c if (isearch_tab(i).eq.1) jj = jj+1
1120 c if (jj.eq.iselect) then
1126 end_select = iran_num(4,nres)
1127 psi=max_phi*gen_psi()
1128 if (psi.eq.0.0D0) then
1132 phi(end_select)=pinorm(phi(end_select)+psi)
1133 if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)')
1134 & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
1135 & phi(end_select)*rad2deg
1136 c if (end_select.gt.3)
1137 c & theta(end_select-1)=gen_theta(itype(end_select-2),
1138 c & phi(end_select-1),phi(end_select))
1139 c if (end_select.lt.nres)
1140 c & theta(end_select)=gen_theta(itype(end_select-1),
1141 c & phi(end_select),phi(end_select+1))
1142 cd print *,'nres=',nres,' end_select=',end_select
1143 cd print *,'theta',end_select-1,theta(end_select-1)
1144 cd print *,'theta',end_select,theta(end_select)
1149 C------------------------------------------------------------------------------
1151 C Select the number of SCs to perturb.
1153 301 nside_move=iran_num(1,MaxSideMove)
1154 c print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
1155 C Select the indices.
1158 111 inds=iran_num(nnt,nct)
1160 if (icount.gt.1000) then
1161 write (iout,'(a)')'Error - cannot select side chains to move.'
1165 if (itype(inds).eq.10) goto 111
1167 if (inds.eq.ind_side(j)) goto 111
1170 if (inds.lt.ind_side(j)) then
1176 112 do j=i,indx+1,-1
1177 ind_side(j)=ind_side(j-1)
1179 113 ind_side(indx)=inds
1181 C Carry out perturbation.
1185 call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
1188 if (isctry.gt.1000) then
1189 write (iout,'(a)') 'Too many errors in SC generation.'
1195 if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)')
1196 & 'Side chain ',restyp(iti),ii,' moved to ',
1197 & alph(ii)*rad2deg,omeg(ii)*rad2deg
1202 C------------------------------------------------------------------------------
1204 400 end_select=iran_num(3,nres)
1205 theta_new=gen_theta(itype(end_select),phi(end_select),
1206 & phi(end_select+1))
1207 if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)')
1208 & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
1209 & ' to ',theta_new*rad2deg
1210 theta(end_select)=theta_new
1214 C------------------------------------------------------------------------------
1215 C Error returned from SelectMove.
1219 C------------------------------------------------------------------------------
1220 subroutine SelectMove(MoveType)
1221 implicit real*8 (a-h,o-z)
1222 include 'DIMENSIONS'
1223 include 'COMMON.MCM'
1224 include 'COMMON.IOUNITS'
1225 what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
1227 if (what_move.ge.sumpro_type(i-1).and.
1228 & what_move.lt.sumpro_type(i)) then
1234 & 'Fatal error in SelectMoveType: cannot select move.'
1235 MoveType=MaxMoveType+1
1238 c----------------------------------------------------------------------------
1239 double precision function gen_psi()
1242 double precision x,ran_number
1243 include 'COMMON.IOUNITS'
1244 include 'COMMON.GEO'
1247 x=ran_number(-pi,pi)
1248 if (dabs(x).gt.angmin) then
1253 write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
1257 c----------------------------------------------------------------------------
1258 subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
1260 implicit real*8 (a-h,o-z)
1261 include 'DIMENSIONS'
1262 include 'COMMON.MCM'
1263 include 'COMMON.IOUNITS'
1264 include 'COMMON.VAR'
1265 include 'COMMON.GEO'
1266 crc include 'COMMON.DEFORM'
1267 double precision ecur,eold,xx,ran_number,bol
1268 double precision xcur(n),xold(n)
1269 double precision ecut1 ,ecut2 ,tola
1270 logical accepted,similar,not_done,enelower
1272 data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
1276 C Set lprn=.true. for debugging.
1279 &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
1283 C Check if the conformation is similar.
1285 reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
1287 write (iout,*) 'Metropolis'
1288 write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
1290 C If energy went down remarkably, we accept the new conformation
1292 cjp if (reldife.lt.ecut1) then
1293 if (difene.lt.ecut1) then
1296 if (lprn) write (iout,'(a)')
1297 & 'Conformation accepted, because energy has lowered remarkably.'
1298 ! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
1299 cjp elseif (reldife.lt.ecut2)
1300 elseif (difene.lt.ecut2)
1302 C Reject the conf. if energy has changed insignificantly and there is not
1303 C much change in conformation.
1305 & write (iout,'(2a)') 'Conformation rejected, because it is',
1306 & ' similar to the preceding one.'
1310 C Else carry out Metropolis test.
1312 xx=ran_number(0.0D0,1.0D0)
1315 & write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
1316 if (xxh.gt.50.0D0) then
1321 if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
1324 if (lprn) write (iout,'(a)')
1325 & 'Conformation accepted, because it passed Metropolis test.'
1328 if (lprn) write (iout,'(a)')
1329 & 'Conformation rejected, because it did not pass Metropolis test.'
1339 c------------------------------------------------------------------------------
1340 integer function conf_comp(x,ene)
1341 implicit real*8 (a-h,o-z)
1342 include 'DIMENSIONS'
1343 include 'COMMON.MCM'
1344 include 'COMMON.VAR'
1345 include 'COMMON.IOUNITS'
1346 include 'COMMON.GEO'
1347 double precision etol , angtol
1348 double precision x(maxvar)
1349 double precision dif_ang,difa
1350 data etol /0.1D0/, angtol /20.0D0/
1352 c write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
1353 if (dabs(ene-esave(ii)).lt.etol) then
1354 difa=dif_ang(nphi,x,varsave(1,ii))
1356 c write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1357 c & rad2deg*varsave(i,ii)
1359 c write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
1360 if (difa.le.angtol) then
1361 if (print_mc.gt.0) then
1362 write (iout,'(a,i5,2(a,1pe15.4))')
1363 & 'Current conformation matches #',ii,
1364 & ' in the store array ene=',ene,' esave=',esave(ii)
1365 c write (*,'(a,i5,a)') 'Current conformation matches #',ii,
1366 c & ' in the store array.'
1367 endif ! print_mc.gt.0
1368 if (print_mc.gt.1) then
1370 write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1371 & rad2deg*varsave(i,ii)
1373 endif ! print_mc.gt.1
1383 C----------------------------------------------------------------------------
1384 double precision function dif_ang(n,x,y)
1387 double precision x(n),y(n)
1388 double precision w,wa,dif,difa
1389 double precision pinorm
1390 include 'COMMON.GEO'
1394 dif=dabs(pinorm(y(i)-x(i)))
1395 if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
1396 w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
1400 dif_ang=rad2deg*dsqrt(difa/wa)
1403 c--------------------------------------------------------------------------
1404 subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
1405 & ecur,xcur,ecache,xcache)
1407 include 'COMMON.GEO'
1408 include 'COMMON.IOUNITS'
1409 integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
1411 double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
1412 cd write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
1413 cd write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
1414 cd write (iout,*) 'Old CACHE array:'
1416 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1417 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1421 do while (i.gt.0 .and. ecur.lt.ecache(i))
1425 cd write (iout,*) 'i=',i,' ncache=',ncache
1426 if (ncache.eq.n2) then
1427 write (iout,*) 'Cache dimension exceeded',ncache,n2
1428 write (iout,*) 'Highest-energy conformation will be removed.'
1432 ecache(ii+1)=ecache(ii)
1433 CachSrc(ii+1)=CachSrc(ii)
1435 xcache(j,ii+1)=xcache(j,ii)
1444 cd write (iout,*) 'New CACHE array:'
1446 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1447 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1451 c--------------------------------------------------------------------------
1452 subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
1455 include 'COMMON.GEO'
1456 include 'COMMON.IOUNITS'
1457 integer n1,n2,ncache,nvar,CachSrc(n2)
1459 double precision ecache(n2),xcache(n1,n2)
1461 cd write (iout,*) 'Enter RM_FROM_CACHE'
1462 cd write (iout,*) 'Old CACHE array:'
1464 cd write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
1465 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
1469 ecache(ii-1)=ecache(ii)
1470 CachSrc(ii-1)=CachSrc(ii)
1472 xcache(j,ii-1)=xcache(j,ii)
1476 cd write (iout,*) 'New CACHE array:'
1478 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1479 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)