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,1),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
262 call gen_rand_conf(nstart_grow,nrestmp,*30)
264 call geom_to_var(nvar,varia)
265 cd write (iout,'(a)') 'New variables:'
266 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
269 call etotal(energia(0))
271 c call enerprint(energia(0))
272 c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
273 if (etot-elowest.gt.overlap_cut) then
274 if(iprint.gt.1.or.etot.lt.1d20)
275 & write (iout,'(a,1pe14.5)')
276 & 'Overlap detected in the current conf.; energy is',etot
280 if (noverlap.gt.maxoverlap) then
281 write (iout,'(a)') 'Too many overlapping confs.'
286 call minimize(etot,varia,iretcode,nfun)
287 cd write (iout,*) 'etot from MINIMIZE:',etot
288 cd write (iout,'(a)') 'Variables after minimization:'
289 cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
291 call etotal(energia(0))
293 neneval=neneval+nfun+2
295 c call enerprint(energia(0))
296 write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,
298 C--------------------------------------------------------------------------
299 C... Do Metropolis test
300 C--------------------------------------------------------------------------
304 if (WhatsUp.eq.0 .and. Kwita.eq.0) then
305 call metropolis(nvar,varia,varold,etot,eold,accepted,
306 & my_conf,EneLower,it)
308 write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
313 if (elowest.gt.etot) elowest=etot
314 moves_acc(MoveType)=moves_acc(MoveType)+1
315 if (MoveType.eq.1) then
316 nbond_acc(nbond)=nbond_acc(nbond)+1
318 C Check against conformation repetitions.
319 irepet=conf_comp(varia,etot)
321 #if defined(AIX) || defined(PGI) || defined(CRAY)
322 open (istat,file=statname,position='append')
324 open (istat,file=statname,access='append')
327 call statprint(nacc,nfun,iretcode,etot,elowest)
329 call var_to_geom(nvar,varia)
331 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),
332 & nsup,przes,obr,non_conv)
334 call contact(.false.,ncont,icont,co)
335 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
336 write (iout,'(a,f8.3,a,f8.3)')
337 & 'RMS deviation from the reference structure:',rms,
338 & ' % of native contacts:',frac*100,' contact order',co
342 write (iout,*) 'Writing new conformation',nout
344 write (istat,'(i5,16(1pe14.5))') nout,
345 & (energia(print_order(i)),i=1,nprint_ene),
349 & write (istat,'(i5,17(1pe14.5))') nout,
350 & (energia(print_order(i)),i=1,nprint_ene),etot
352 if (print_stat) close(istat)
353 C Print internal coordinates.
354 if (print_int) call briefout(nout,etot)
355 C Accumulate the newly accepted conf in the coord1 array, if it is different
356 C from all confs that are already there.
357 call compare_s1(n_thr,max_thread2,etot,varia,ii,
358 & enetb1,coord1,rms_deform,.true.,iprint)
359 write (iout,*) 'After compare_ss: n_thr=',n_thr
360 if (ii.eq.1 .or. ii.eq.3) then
361 write (iout,'(8f10.4)')
362 & (rad2deg*coord1(i,n_thr),i=1,nvar)
365 write (iout,*) 'Conformation from cache, not written.'
368 if (nrepm.gt.maxrepm) then
369 write (iout,'(a)') 'Too many conformation repetitions.'
372 C Store the accepted conf. and its energy.
377 if (irepet.eq.0) call zapis(varia,etot)
378 C Lower the temperature, if necessary.
389 C Check for time limit.
390 if (ovrtim()) WhatsUp=-1
391 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
400 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
401 call statprint(nacc,nfun,iretcode,etot,elowest)
403 & 'Statistics of multiple-bond motions. Total motions:'
404 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
405 write (iout,'(a)') 'Accepted motions:'
406 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
408 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
414 c------------------------------------------------------------------------------
415 subroutine do_mcm(i_orig)
416 C Monte-Carlo-with-Minimization calculations - parallel code.
417 implicit real*8 (a-h,o-z)
420 include 'COMMON.IOUNITS'
422 include 'COMMON.CHAIN'
424 include 'COMMON.CONTACTS'
425 include 'COMMON.CONTROL'
427 include 'COMMON.INTERACT'
428 include 'COMMON.INFO'
429 include 'COMMON.CACHE'
430 crc include 'COMMON.DEFORM'
431 crc include 'COMMON.DEFORM1'
432 crc include 'COMMON.DEFORM2'
433 include 'COMMON.MINIM'
434 include 'COMMON.NAMES'
435 logical accepted,over,ovrtim,error,lprint,not_done,similar,
436 & enelower,non_conv,flag,finish
437 integer MoveType,nbond,conf_comp
438 double precision varia(maxvar),varold(maxvar),elowest,eold,
439 & x1(maxvar), varold1(maxvar), przes(3),obr(3,3)
440 integer iparentx(max_threadss2)
441 integer iparentx1(max_threadss2)
442 integer imtasks(150),imtasks_n
443 double precision energia(0:n_ene)
445 print *,'Master entered DO_MCM'
453 C---------------------------------------------------------------------------
454 C Initialize counters.
455 C---------------------------------------------------------------------------
456 C Total number of generated confs.
458 C Total number of moves. In general this won`t be equal to the number of
459 C attempted moves, because we may want to reject some "bad" confs just by
462 C Total number of temperature jumps.
464 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
470 C Initialize total and accepted number of moves of various kind.
475 C Total number of energy evaluations.
479 c write (iout,*) 'RanFract=',RanFract
482 c----------------------------------------------------------------------------
483 C Compute and print initial energies.
484 c----------------------------------------------------------------------------
486 write (iout,'(/80(1h*)/a)') 'Initial energies:'
489 call etotal(energia(0))
491 call enerprint(energia(0))
492 C Request energy computation from slave processors.
493 call geom_to_var(nvar,varia)
494 ! write (iout,*) 'The VARIA array'
495 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
496 call minimize(etot,varia,iretcode,nfun)
497 call var_to_geom(nvar,varia)
499 write (iout,*) 'etot from MINIMIZE:',etot
500 ! write (iout,*) 'Tha VARIA array'
501 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
504 if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))')
505 & 'Enter Monte Carlo procedure.'
506 if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
512 call zapis(varia,etot)
514 call var_to_geom(nvar,varia)
516 call etotal(energia(0))
517 if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
519 nacc=0 ! total # of accepted confs of the current processor.
520 nacc_tot=0 ! total # of accepted confs of all processors.
522 C----------------------------------------------------------------------------
524 c----------------------------------------------------------------------------
527 LOOP1:do while (not_done)
529 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
530 & 'Beginning iteration #',it
531 C Initialize local counter.
532 ntrial=0 ! # of generated non-overlapping confs.
533 noverlap=0 ! # of overlapping confs.
535 LOOP2:do while (.not. accepted)
537 LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
539 if(imtasks(i).eq.0) then
544 C Retrieve the angles of previously accepted conformation
548 call var_to_geom(nvar,varia)
551 C Heat up the system, if necessary.
553 C If temperature cannot be further increased, stop.
559 c write (iout,'(a)') 'Old variables:'
560 c write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
561 C Decide whether to generate a random conformation or perturb the old one
562 RandOrPert=ran_number(0.0D0,1.0D0)
563 if (RandOrPert.gt.RanFract) then
565 & write (iout,'(a)') 'Perturbation-generated conformation.'
566 call perturb(error,lprint,MoveType,nbond,1.0D0)
567 c print *,'after perturb',error,finish
568 if (error) finish = .true.
569 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
570 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
571 & MoveType,' returned from PERTURB.'
573 write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',
574 & MoveType,' returned from PERTURB.'
580 nstart_grow=iran_num(3,nres)
582 & write (iout,'(2a,i3)') 'Random-generated conformation',
583 & ' - chain regrown from residue',nstart_grow
585 call gen_rand_conf(nstart_grow,nrestmp,*30)
587 call geom_to_var(nvar,varia)
589 c print *,'finish=',finish
590 if (etot-elowest.gt.overlap_cut) then
591 if (print_mc.gt.1) write (iout,'(a,1pe14.5)')
592 & 'Overlap detected in the current conf.; energy is',etot
593 if(iprint.gt.1.or.etot.lt.1d19) print *,
594 & 'Overlap detected in the current conf.; energy is',etot
598 if (noverlap.gt.maxoverlap) then
599 write (iout,*) 'Too many overlapping confs.',
600 & ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
603 else if (.not. finish) then
604 C Distribute tasks to processors
605 c print *,'Master sending order'
606 call MPI_SEND(12, 1, MPI_INTEGER, is, tag,
608 c write (iout,*) '12: tag=',tag
609 c print *,'Master sent order to processor',is
610 call MPI_SEND(it, 1, MPI_INTEGER, is, tag,
612 c write (iout,*) 'it: tag=',tag
613 call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,
615 c write (iout,*) 'eold: tag=',tag
616 call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION,
619 c write (iout,*) 'varia: tag=',tag
620 call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION,
623 c write (iout,*) 'varold: tag=',tag
630 imtasks_n=imtasks_n+1
636 LOOP_RECV:do while(.not.flag)
638 call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
640 call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,
641 & CG_COMM, status, ierr)
642 call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,
643 & CG_COMM, status, ierr)
644 call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,
645 & CG_COMM, status, ierr)
646 call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,
647 & CG_COMM, status, ierr)
648 call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is,
649 & tag, CG_COMM, status, ierr)
650 call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,
651 & CG_COMM, status, ierr)
652 call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,
653 & CG_COMM, status, ierr)
654 call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,
655 & CG_COMM, status, ierr)
656 i_grnum_d=i_grnum_d+ii_grnum_d
657 i_ennum_d=i_ennum_d+ii_ennum_d
658 neneval = neneval+ii_ennum_d
659 i_hesnum_d=i_hesnum_d+ii_hesnum_d
660 i_minimiz=i_minimiz+1
662 imtasks_n=imtasks_n-1
668 if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)')
669 & 'From Worker #',is,' iitt',iitt,
670 & ' Conformation:',ngen,' energy:',etot
671 C--------------------------------------------------------------------------
672 C... Do Metropolis test
673 C--------------------------------------------------------------------------
674 call metropolis(nvar,varia,varold1,etot,eold1,accepted,
676 if(iitt.ne.it.and..not.similar) then
677 call metropolis(nvar,varia,varold,etot,eold,accepted,
681 if(etot.lt.eneglobal)eneglobal=etot
682 c if(mod(it,100).eq.0)
683 write(iout,*)'CHUJOJEB ',neneval,eneglobal
685 C Write the accepted conformation.
688 call var_to_geom(nvar,varia)
691 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
692 & nsup,przes,obr,non_conv)
694 call contact(.false.,ncont,icont,co)
695 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
696 write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
697 & 'RMS deviation from the reference structure:',rms,
698 & ' % of native contacts:',frac*100,' contact order:',co
701 & write (iout,*) 'Writing new conformation',nout
703 call var_to_geom(nvar,varia)
704 #if defined(AIX) || defined(PGI) || defined(CRAY)
705 open (istat,file=statname,position='append')
707 open (istat,file=statname,access='append')
710 write (istat,'(i5,16(1pe14.5))') nout,
711 & (energia(print_order(i)),i=1,nprint_ene),
714 write (istat,'(i5,16(1pe14.5))') nout,
715 & (energia(print_order(i)),i=1,nprint_ene),etot
719 C Print internal coordinates.
720 if (print_int) call briefout(nout,etot)
723 if (elowest.gt.etot) elowest=etot
724 moves_acc(MoveType)=moves_acc(MoveType)+1
725 if (MoveType.eq.1) then
726 nbond_acc(nbond)=nbond_acc(nbond)+1
728 C Check against conformation repetitions.
729 irepet=conf_comp(varia,etot)
730 if (nrepm.gt.maxrepm) then
732 & write (iout,'(a)') 'Too many conformation repetitions.'
735 C Store the accepted conf. and its energy.
740 if (irepet.eq.0) call zapis(varia,etot)
741 C Lower the temperature, if necessary.
747 if(finish.and.imtasks_n.eq.0)exit LOOP2
748 enddo LOOP2 ! accepted
749 C Check for time limit.
750 not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
751 if(.not.not_done .or. finish) then
752 if(imtasks_n.gt.0) then
759 enddo LOOP1 ! not_done
761 if (print_mc.gt.0) then
762 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
763 call statprint(nacc,nfun,iretcode,etot,elowest)
765 & 'Statistics of multiple-bond motions. Total motions:'
766 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
767 write (iout,'(a)') 'Accepted motions:'
768 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
770 &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
778 call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
783 c------------------------------------------------------------------------------
784 subroutine execute_slave(nodeinfo,iprint)
785 implicit real*8 (a-h,o-z)
788 include 'COMMON.TIME1'
789 include 'COMMON.IOUNITS'
790 crc include 'COMMON.DEFORM'
791 crc include 'COMMON.DEFORM1'
792 crc include 'COMMON.DEFORM2'
793 include 'COMMON.LOCAL'
795 include 'COMMON.INFO'
796 include 'COMMON.MINIM'
797 character*10 nodeinfo
798 double precision x(maxvar),x1(maxvar)
800 c print *,'Processor:',MyID,' Entering execute_slave'
802 c call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
805 1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
806 & CG_COMM, status, ierr)
807 c write(iout,*)'12: tag=',tag
808 if(iprint.ge.2)print *, MyID,' recv order ',i_switch
809 if (i_switch.eq.12) then
813 call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
814 & CG_COMM, status, ierr)
815 c write(iout,*)'12: tag=',tag
816 call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
817 & CG_COMM, status, ierr)
818 c write(iout,*)'ener: tag=',tag
819 call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
820 & CG_COMM, status, ierr)
821 c write(iout,*)'x: tag=',tag
822 call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
823 & CG_COMM, status, ierr)
824 c write(iout,*)'x1: tag=',tag
830 c print *,'calling minimize'
831 call minimize(energyx,x,iretcode,nfun)
833 & write(iout,100)'minimized energy = ',energyx,
834 & ' # funeval:',nfun,' iret ',iretcode
835 write(*,100)'minimized energy = ',energyx,
836 & ' # funeval:',nfun,' iret ',iretcode
837 100 format(a20,f10.5,a12,i5,a6,i2)
838 if(iretcode.eq.10) then
841 & write(iout,*)' ... not converged - trying again ',iminrep
842 call minimize(energyx,x,iretcode,nfun)
844 & write(iout,*)'minimized energy = ',energyx,
845 & ' # funeval:',nfun,' iret ',iretcode
846 if(iretcode.ne.10)go to 812
848 if(iretcode.eq.10) then
850 & write(iout,*)' ... not converged again - giving up'
855 c print *,'Sending results'
856 call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
858 call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
860 call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
862 call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
864 call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
866 call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
868 call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
870 call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
872 c print *,'End sending'
879 c------------------------------------------------------------------------------
880 subroutine statprint(it,nfun,iretcode,etot,elowest)
881 implicit real*8 (a-h,o-z)
883 include 'COMMON.IOUNITS'
884 include 'COMMON.CONTROL'
888 & '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
889 & 'Finished iteration #',it,' energy is',etot,
890 & ' lowest energy:',elowest,
891 & 'SUMSL return code:',iretcode,
892 & ' # of energy evaluations:',neneval,
893 & '# of temperature jumps:',ntherm,
894 & ' # of minima repetitions:',nrepm
896 write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
897 & 'Finished iteration #',it,' energy is',etot,
898 & ' lowest energy:',elowest
901 & 'Kind of move ',' total',' accepted',
903 write (iout,'(58(1h-))')
905 if (moves(i).eq.0) then
908 fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
910 write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
913 write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,
914 & dfloat(nacc_tot)/dfloat(nmove)
915 write (iout,'(58(1h-))')
916 write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
919 c------------------------------------------------------------------------------
920 subroutine heat(over)
921 implicit real*8 (a-h,o-z)
924 include 'COMMON.IOUNITS'
926 C Check if there`s a need to increase temperature.
927 if (ntrial.gt.maxtrial) then
928 if (NstepH.gt.0) then
929 if (dabs(Tcur-TMax).lt.1.0D-7) then
931 & write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))')
932 & 'Upper limit of temperature reached. Terminating.'
937 if (Tcur.gt.Tmax) Tcur=Tmax
938 betbol=1.0D0/(Rbol*Tcur)
940 & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
941 & 'System heated up to ',Tcur,' K; BetBol:',betbol
949 & 'Maximum number of trials in a single MCM iteration exceeded.'
958 c------------------------------------------------------------------------------
960 implicit real*8 (a-h,o-z)
963 include 'COMMON.IOUNITS'
964 if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
966 if (Tcur.lt.Tmin) Tcur=Tmin
967 betbol=1.0D0/(Rbol*Tcur)
969 & write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
970 & 'System cooled down up to ',Tcur,' K; BetBol:',betbol
974 C---------------------------------------------------------------------------
975 subroutine zapis(varia,etot)
976 implicit real*8 (a-h,o-z)
980 include 'COMMON.INFO'
985 include 'COMMON.IOUNITS'
986 integer itemp(maxsave)
987 double precision varia(maxvar)
991 write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
992 & ' MaxSave=',MaxSave
993 write (iout,'(a)') 'Current energy and conformation:'
994 write (iout,'(1pe14.5)') etot
995 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
997 C Shift the contents of the esave and varsave arrays if filled up.
998 call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
999 & etot,varia,esave,varsave)
1001 write (iout,'(a)') 'Energies and the VarSave array.'
1003 write (iout,'(i5,1pe14.5)') i,esave(i)
1004 write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
1009 C---------------------------------------------------------------------------
1010 subroutine perturb(error,lprint,MoveType,nbond,max_phi)
1011 implicit real*8 (a-h,o-z)
1012 include 'DIMENSIONS'
1013 parameter (MMaxSideMove=100)
1014 include 'COMMON.MCM'
1015 include 'COMMON.CHAIN'
1016 include 'COMMON.INTERACT'
1017 include 'COMMON.VAR'
1018 include 'COMMON.GEO'
1019 include 'COMMON.NAMES'
1020 include 'COMMON.IOUNITS'
1021 crc include 'COMMON.DEFORM1'
1022 logical error,lprint,fail
1023 integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
1024 double precision max_phi
1025 double precision psi,gen_psi
1032 C Perturb the conformation according to a randomly selected move.
1033 call SelectMove(MoveType)
1034 c write (iout,*) 'MoveType=',MoveType
1036 goto (100,200,300,400,500) MoveType
1037 C------------------------------------------------------------------------------
1038 C Backbone N-bond move.
1039 C Select the number of bonds (length of the segment to perturb).
1041 if (itrial.gt.1000) then
1042 write (iout,'(a)') 'Too many attempts at multiple-bond move.'
1046 bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
1047 c print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
1048 c & ' Bond_prob=',Bond_Prob
1050 c print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
1051 if (bond_prob.ge.sumpro_bond(i) .and.
1052 & bond_prob.le.sumpro_bond(i+1)) then
1057 write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
1058 & ' to move out of range.'
1062 if (nwindow.gt.0) then
1063 C Select the first residue to perturb
1064 iwindow=iran_num(1,nwindow)
1065 print *,'iwindow=',iwindow
1067 do while (winlen(iwindow).lt.nbond)
1068 iwindow=iran_num(1,nwindow)
1070 if (iiwin.gt.1000) then
1071 write (iout,'(a)') 'Cannot select moveable residues.'
1076 nstart=iran_num(winstart(iwindow),winend(iwindow))
1078 nstart = iran_num(koniecl+2,nres-nbond-koniecl)
1079 cd print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
1080 cd & ' nstart=',nstart
1083 if (psi.eq.0.0) then
1087 if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
1088 & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
1089 cd print *,'nstart=',nstart
1090 call bond_move(nbond,nstart,psi,.false.,error)
1093 & 'Could not define reference system in bond_move, ',
1094 & 'choosing ahother segment.'
1098 nbond_move(nbond)=nbond_move(nbond)+1
1102 C------------------------------------------------------------------------------
1103 C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
1107 c end_select=iran_num(1,2*koniecl)
1108 c if (end_select.gt.koniecl) then
1109 c end_select=nphi-(end_select-koniecl)
1111 c end_select=koniecl+3
1113 c if (nwindow.gt.0) then
1114 c iwin=iran_num(1,nwindow)
1115 c i1=max0(4,winstart(iwin))
1116 c i2=min0(winend(imin)+2,nres)
1117 c end_select=iran_num(i1,i2)
1119 c iselect = iran_num(1,nmov_var)
1122 c if (isearch_tab(i).eq.1) jj = jj+1
1123 c if (jj.eq.iselect) then
1129 end_select = iran_num(4,nres)
1130 psi=max_phi*gen_psi()
1131 if (psi.eq.0.0D0) then
1135 phi(end_select)=pinorm(phi(end_select)+psi)
1136 if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)')
1137 & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
1138 & phi(end_select)*rad2deg
1139 c if (end_select.gt.3)
1140 c & theta(end_select-1)=gen_theta(itype(end_select-2),
1141 c & phi(end_select-1),phi(end_select))
1142 c if (end_select.lt.nres)
1143 c & theta(end_select)=gen_theta(itype(end_select-1),
1144 c & phi(end_select),phi(end_select+1))
1145 cd print *,'nres=',nres,' end_select=',end_select
1146 cd print *,'theta',end_select-1,theta(end_select-1)
1147 cd print *,'theta',end_select,theta(end_select)
1152 C------------------------------------------------------------------------------
1154 C Select the number of SCs to perturb.
1156 301 nside_move=iran_num(1,MaxSideMove)
1157 c print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
1158 C Select the indices.
1161 111 inds=iran_num(nnt,nct)
1163 if (icount.gt.1000) then
1164 write (iout,'(a)')'Error - cannot select side chains to move.'
1168 if (itype(inds).eq.10) goto 111
1170 if (inds.eq.ind_side(j)) goto 111
1173 if (inds.lt.ind_side(j)) then
1179 112 do j=i,indx+1,-1
1180 ind_side(j)=ind_side(j-1)
1182 113 ind_side(indx)=inds
1184 C Carry out perturbation.
1188 call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
1191 if (isctry.gt.1000) then
1192 write (iout,'(a)') 'Too many errors in SC generation.'
1198 if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)')
1199 & 'Side chain ',restyp(iti),ii,' moved to ',
1200 & alph(ii)*rad2deg,omeg(ii)*rad2deg
1205 C------------------------------------------------------------------------------
1207 400 end_select=iran_num(3,nres)
1208 theta_new=gen_theta(itype(end_select),phi(end_select),
1209 & phi(end_select+1))
1210 if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)')
1211 & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
1212 & ' to ',theta_new*rad2deg
1213 theta(end_select)=theta_new
1217 C------------------------------------------------------------------------------
1218 C Error returned from SelectMove.
1222 C------------------------------------------------------------------------------
1223 subroutine SelectMove(MoveType)
1224 implicit real*8 (a-h,o-z)
1225 include 'DIMENSIONS'
1226 include 'COMMON.MCM'
1227 include 'COMMON.IOUNITS'
1228 what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
1230 if (what_move.ge.sumpro_type(i-1).and.
1231 & what_move.lt.sumpro_type(i)) then
1237 & 'Fatal error in SelectMoveType: cannot select move.'
1238 MoveType=MaxMoveType+1
1241 c----------------------------------------------------------------------------
1242 double precision function gen_psi()
1245 double precision x,ran_number
1246 include 'COMMON.IOUNITS'
1247 include 'COMMON.GEO'
1250 x=ran_number(-pi,pi)
1251 if (dabs(x).gt.angmin) then
1256 write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
1260 c----------------------------------------------------------------------------
1261 subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
1263 implicit real*8 (a-h,o-z)
1264 include 'DIMENSIONS'
1265 include 'COMMON.MCM'
1266 include 'COMMON.IOUNITS'
1267 include 'COMMON.VAR'
1268 include 'COMMON.GEO'
1269 crc include 'COMMON.DEFORM'
1270 double precision ecur,eold,xx,ran_number,bol
1271 double precision xcur(n),xold(n)
1272 double precision ecut1 ,ecut2 ,tola
1273 logical accepted,similar,not_done,enelower
1275 data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
1279 C Set lprn=.true. for debugging.
1282 &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
1286 C Check if the conformation is similar.
1288 reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
1290 write (iout,*) 'Metropolis'
1291 write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
1293 C If energy went down remarkably, we accept the new conformation
1295 cjp if (reldife.lt.ecut1) then
1296 if (difene.lt.ecut1) then
1299 if (lprn) write (iout,'(a)')
1300 & 'Conformation accepted, because energy has lowered remarkably.'
1301 ! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
1302 cjp elseif (reldife.lt.ecut2)
1303 elseif (difene.lt.ecut2)
1305 C Reject the conf. if energy has changed insignificantly and there is not
1306 C much change in conformation.
1308 & write (iout,'(2a)') 'Conformation rejected, because it is',
1309 & ' similar to the preceding one.'
1313 C Else carry out Metropolis test.
1315 xx=ran_number(0.0D0,1.0D0)
1318 & write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
1319 if (xxh.gt.50.0D0) then
1324 if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
1327 if (lprn) write (iout,'(a)')
1328 & 'Conformation accepted, because it passed Metropolis test.'
1331 if (lprn) write (iout,'(a)')
1332 & 'Conformation rejected, because it did not pass Metropolis test.'
1342 c------------------------------------------------------------------------------
1343 integer function conf_comp(x,ene)
1344 implicit real*8 (a-h,o-z)
1345 include 'DIMENSIONS'
1346 include 'COMMON.MCM'
1347 include 'COMMON.VAR'
1348 include 'COMMON.IOUNITS'
1349 include 'COMMON.GEO'
1350 double precision etol , angtol
1351 double precision x(maxvar)
1352 double precision dif_ang,difa
1353 data etol /0.1D0/, angtol /20.0D0/
1355 c write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
1356 if (dabs(ene-esave(ii)).lt.etol) then
1357 difa=dif_ang(nphi,x,varsave(1,ii))
1359 c write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1360 c & rad2deg*varsave(i,ii)
1362 c write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
1363 if (difa.le.angtol) then
1364 if (print_mc.gt.0) then
1365 write (iout,'(a,i5,2(a,1pe15.4))')
1366 & 'Current conformation matches #',ii,
1367 & ' in the store array ene=',ene,' esave=',esave(ii)
1368 c write (*,'(a,i5,a)') 'Current conformation matches #',ii,
1369 c & ' in the store array.'
1370 endif ! print_mc.gt.0
1371 if (print_mc.gt.1) then
1373 write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1374 & rad2deg*varsave(i,ii)
1376 endif ! print_mc.gt.1
1386 C----------------------------------------------------------------------------
1387 double precision function dif_ang(n,x,y)
1390 double precision x(n),y(n)
1391 double precision w,wa,dif,difa
1392 double precision pinorm
1393 include 'COMMON.GEO'
1397 dif=dabs(pinorm(y(i)-x(i)))
1398 if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
1399 w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
1403 dif_ang=rad2deg*dsqrt(difa/wa)
1406 c--------------------------------------------------------------------------
1407 subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
1408 & ecur,xcur,ecache,xcache)
1410 include 'COMMON.GEO'
1411 include 'COMMON.IOUNITS'
1412 integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
1414 double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
1415 cd write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
1416 cd write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
1417 cd write (iout,*) 'Old CACHE array:'
1419 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1420 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1424 do while (i.gt.0 .and. ecur.lt.ecache(i))
1428 cd write (iout,*) 'i=',i,' ncache=',ncache
1429 if (ncache.eq.n2) then
1430 write (iout,*) 'Cache dimension exceeded',ncache,n2
1431 write (iout,*) 'Highest-energy conformation will be removed.'
1435 ecache(ii+1)=ecache(ii)
1436 CachSrc(ii+1)=CachSrc(ii)
1438 xcache(j,ii+1)=xcache(j,ii)
1447 cd write (iout,*) 'New CACHE array:'
1449 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1450 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1454 c--------------------------------------------------------------------------
1455 subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
1458 include 'COMMON.GEO'
1459 include 'COMMON.IOUNITS'
1460 integer n1,n2,ncache,nvar,CachSrc(n2)
1462 double precision ecache(n2),xcache(n1,n2)
1464 cd write (iout,*) 'Enter RM_FROM_CACHE'
1465 cd write (iout,*) 'Old CACHE array:'
1467 cd write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
1468 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
1472 ecache(ii-1)=ecache(ii)
1473 CachSrc(ii-1)=CachSrc(ii)
1475 xcache(j,ii-1)=xcache(j,ii)
1479 cd write (iout,*) 'New CACHE array:'
1481 cd write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1482 cd write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)