wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / mcm.F
1       subroutine mcm_setup
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.IOUNITS'
5       include 'COMMON.MCM'
6       include 'COMMON.CONTROL'
7       include 'COMMON.INTERACT'
8       include 'COMMON.NAMES'
9       include 'COMMON.CHAIN'
10       include 'COMMON.VAR'
11 C
12 C Set up variables used in MC/MCM.
13 C    
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:'
22         do i=1,nwindow
23           i1=winstart(i)
24           i2=winend(i)
25           it1=itype(i1)
26           it2=itype(i2)
27           write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,
28      &                        ' length',winlen(i)
29         enddo
30       endif
31 C Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
32       RBol=1.9858D-3
33 C Number of "end bonds".
34       koniecl=0
35 c     koniecl=nphi
36       print *,'koniecl=',koniecl
37       write (iout,'(a)') 'Probabilities of move types:'
38       write (*,'(a)') 'Probabilities of move types:'
39       do i=1,MaxMoveType
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)
44       enddo
45       write (iout,*) 
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
49         maxwinlen=winlen(1)
50         do i=2,nwindow
51           if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
52         enddo
53         nbm=min0(maxwinlen,6)
54         write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
55       else
56         nbm=min0(6,nres-2)
57       endif
58       sumpro_bond(0)=0.0D0
59       sumpro_bond(1)=0.0D0 
60       do i=2,nbm
61         sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
62       enddo
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
69       ngly=0
70       do i=nnt,nct
71         if (itype(i).eq.10) ngly=ngly+1
72       enddo
73       mmm=nct-nnt-ngly+1
74       if (mmm.gt.0) then
75         MaxSideMove=min0((nct-nnt+1)/2,mmm)
76       endif
77 c     print *,'MaxSideMove=',MaxSideMove
78 C Max. number of generated confs (not used at present).
79       maxgen=10000
80 C Set initial temperature
81       Tcur=Tmin
82       betbol=1.0D0/(Rbol*Tcur)
83       write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,
84      &    ' BetBol:',betbol
85       write (iout,*) 'RanFract=',ranfract
86       return
87       end
88 c------------------------------------------------------------------------------
89 #ifndef MPI
90       subroutine do_mcm(i_orig)
91 C Monte-Carlo-with-Minimization calculations - serial code.
92       implicit real*8 (a-h,o-z)
93       include 'DIMENSIONS'
94       include 'COMMON.IOUNITS'
95       include 'COMMON.GEO'
96       include 'COMMON.CHAIN'
97       include 'COMMON.MCM'
98       include 'COMMON.CONTACTS'
99       include 'COMMON.CONTROL'
100       include 'COMMON.VAR'
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,
107      &        enelower,non_conv
108       integer MoveType,nbond,conf_comp
109       integer ifeed(max_cache)
110       double precision varia(maxvar),varold(maxvar),elowest,eold,
111      & przes(3),obr(3,3)
112       double precision energia(0:n_ene)
113       double precision coord1(maxres,3)
114
115
116 C---------------------------------------------------------------------------
117 C Initialize counters.
118 C---------------------------------------------------------------------------
119 C Total number of generated confs.
120       ngen=0
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
123 C overlap check.
124       nmove=0
125 C Total number of temperature jumps.
126       ntherm=0
127 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
128 C motions.
129       ncache=0
130       do i=1,nres
131         nbond_move(i)=0
132       enddo
133 C Initialize total and accepted number of moves of various kind.
134       do i=0,MaxMoveType
135         moves(i)=0
136         moves_acc(i)=0
137       enddo
138 C Total number of energy evaluations.
139       neneval=0
140       nfun=0
141       nsave=0
142
143       write (iout,*) 'RanFract=',RanFract
144
145       WhatsUp=0
146       Kwita=0
147
148 c----------------------------------------------------------------------------
149 C Compute and print initial energies.
150 c----------------------------------------------------------------------------
151       call intout
152       write (iout,'(/80(1h*)/a)') 'Initial energies:'
153       call chainbuild
154       nf=0
155
156       call etotal(energia(0))
157       etot = energia(0)
158 C Minimize the energy of the first conformation.
159       if (minim) then
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)
165         call chainbuild
166         write (iout,*) 'etot from MINIMIZE:',etot
167 !       write (iout,*) 'Tha VARIA array'       
168 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
169
170         call etotal(energia(0))
171         etot=energia(0)
172         call enerprint(energia(0))
173       endif
174       if (refstr) then
175         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,
176      &             obr,non_conv)
177         rms=dsqrt(rms)
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
183         if (print_stat)
184      &  write (istat,'(i5,17(1pe14.5))') 0,
185      &   (energia(print_order(i)),i=1,nprint_ene),
186      &   etot,rms,frac,co
187       else
188         if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,
189      &   (energia(print_order(i)),i=1,nprint_ene),etot
190       endif
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.'
195       if (print_int) then
196         close(igeom)
197         call briefout(0,etot)
198       endif
199       eold=etot
200       do i=1,nvar
201         varold(i)=varia(i)
202       enddo
203       elowest=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.
207
208       not_done = (iretcode.ne.11)
209
210 C----------------------------------------------------------------------------
211 C Main loop.
212 c----------------------------------------------------------------------------
213       it=0
214       nout=0
215       do while (not_done)
216         it=it+1
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.
221         accepted=.false.
222         do while (.not. accepted)
223
224 C Retrieve the angles of previously accepted conformation
225           noverlap=0 ! # of overlapping confs.
226           do j=1,nvar
227             varia(j)=varold(j)
228           enddo
229           call var_to_geom(nvar,varia)
230 C Rebuild the chain.
231           call chainbuild
232 C Heat up the system, if necessary.
233           call heat(over)
234 C If temperature cannot be further increased, stop.
235           if (over) goto 20
236           MoveType=0
237           nbond=0
238           lprint=.true.
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
244             if (print_mc.gt.0)
245      &        write (iout,'(a)') 'Perturbation-generated conformation.'
246             call perturb(error,lprint,MoveType,nbond,1.0D0)
247             if (error) goto 20
248             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
249               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
250      &           MoveType,' returned from PERTURB.'
251               goto 20
252             endif
253             call chainbuild
254           else
255             MoveType=0
256             moves(0)=moves(0)+1
257             nstart_grow=iran_num(3,nres)
258             if (print_mc.gt.0)
259      &        write (iout,'(2a,i3)') 'Random-generated conformation',
260      &        ' - chain regrown from residue',nstart_grow
261             call gen_rand_conf(nstart_grow,*30)
262           endif
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)
266           ngen=ngen+1
267
268           call etotal(energia(0))
269           etot=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
276             neneval=neneval+1 
277             accepted=.false.
278             noverlap=noverlap+1
279             if (noverlap.gt.maxoverlap) then
280               write (iout,'(a)') 'Too many overlapping confs.'
281               goto 20
282             endif
283           else
284             if (minim) then
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)
289
290               call etotal(energia(0))
291               etot = energia(0)
292               neneval=neneval+nfun+2
293             endif
294 c           call enerprint(energia(0))
295             write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,
296      &      ' energy:',etot
297 C--------------------------------------------------------------------------
298 C... Do Metropolis test
299 C--------------------------------------------------------------------------
300             accepted=.false.
301             my_conf=.false.
302
303             if (WhatsUp.eq.0 .and. Kwita.eq.0) then
304               call metropolis(nvar,varia,varold,etot,eold,accepted,
305      &                      my_conf,EneLower,it)
306             endif
307             write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
308             if (accepted) then
309
310               nacc=nacc+1
311               nacc_tot=nacc_tot+1
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
316               endif
317 C Check against conformation repetitions.
318               irepet=conf_comp(varia,etot)
319               if (print_stat) then
320 #if defined(AIX) || defined(PGI)
321               open (istat,file=statname,position='append')
322 #else
323                open (istat,file=statname,access='append')
324 #endif
325               endif
326               call statprint(nacc,nfun,iretcode,etot,elowest)
327               if (refstr) then
328                 call var_to_geom(nvar,varia)
329                 call chainbuild
330                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),
331      &                    nsup,przes,obr,non_conv)
332                 rms=dsqrt(rms)
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
338               endif ! refstr
339               if (My_Conf) then
340                 nout=nout+1
341                 write (iout,*) 'Writing new conformation',nout
342                 if (refstr) then
343                   write (istat,'(i5,16(1pe14.5))') nout,
344      &             (energia(print_order(i)),i=1,nprint_ene),
345      &             etot,rms,frac
346                 else
347                   if (print_stat)
348      &             write (istat,'(i5,17(1pe14.5))') nout,
349      &              (energia(print_order(i)),i=1,nprint_ene),etot
350                 endif ! refstr
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)
362                 endif
363               else
364                 write (iout,*) 'Conformation from cache, not written.'
365               endif ! My_Conf 
366
367               if (nrepm.gt.maxrepm) then
368                 write (iout,'(a)') 'Too many conformation repetitions.'
369                 goto 20
370               endif
371 C Store the accepted conf. and its energy.
372               eold=etot
373               do i=1,nvar
374                 varold(i)=varia(i)
375               enddo
376               if (irepet.eq.0) call zapis(varia,etot)
377 C Lower the temperature, if necessary.
378               call cool
379
380             else
381
382               ntrial=ntrial+1
383             endif ! accepted
384           endif ! overlap
385
386    30     continue
387         enddo ! accepted
388 C Check for time limit.
389         if (ovrtim()) WhatsUp=-1
390         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
391      &       .and. (Kwita.eq.0)
392
393       enddo ! not_done
394       goto 21
395    20 WhatsUp=-3
396
397    21 continue
398       runtime=tcpu()
399       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
400       call statprint(nacc,nfun,iretcode,etot,elowest)
401       write (iout,'(a)') 
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)
406       if (it.ge.maxacc)
407      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
408
409       return
410       end
411 #endif
412 #ifdef MPI
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)
417       include 'DIMENSIONS'
418       include 'mpif.h'
419       include 'COMMON.IOUNITS'
420       include 'COMMON.GEO'
421       include 'COMMON.CHAIN'
422       include 'COMMON.MCM'
423       include 'COMMON.CONTACTS'
424       include 'COMMON.CONTROL'
425       include 'COMMON.VAR'
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)
443
444       print *,'Master entered DO_MCM'
445       nodenum = nprocs
446       
447       finish=.false.
448       imtasks_n=0
449       do i=1,nodenum-1
450        imtasks(i)=0
451       enddo
452 C---------------------------------------------------------------------------
453 C Initialize counters.
454 C---------------------------------------------------------------------------
455 C Total number of generated confs.
456       ngen=0
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
459 C overlap check.
460       nmove=0
461 C Total number of temperature jumps.
462       ntherm=0
463 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
464 C motions.
465       ncache=0
466       do i=1,nres
467         nbond_move(i)=0
468       enddo
469 C Initialize total and accepted number of moves of various kind.
470       do i=0,MaxMoveType
471         moves(i)=0
472         moves_acc(i)=0
473       enddo
474 C Total number of energy evaluations.
475       neneval=0
476       nfun=0
477       nsave=0
478 c      write (iout,*) 'RanFract=',RanFract
479       WhatsUp=0
480       Kwita=0
481 c----------------------------------------------------------------------------
482 C Compute and print initial energies.
483 c----------------------------------------------------------------------------
484       call intout
485       write (iout,'(/80(1h*)/a)') 'Initial energies:'
486       call chainbuild
487       nf=0
488       call etotal(energia(0))
489       etot = 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)
497       call chainbuild
498       write (iout,*) 'etot from MINIMIZE:',etot
499 !     write (iout,*) 'Tha VARIA array'
500 !     write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
501       neneval=0
502       eneglobal=1.0d99
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
506       eold=etot
507       do i=1,nvar
508         varold(i)=varia(i)
509       enddo
510       elowest=etot
511       call zapis(varia,etot)
512 c diagnostics
513       call var_to_geom(nvar,varia)
514       call chainbuild
515       call etotal(energia(0))
516       if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
517 c end diagnostics
518       nacc=0         ! total # of accepted confs of the current processor.
519       nacc_tot=0     ! total # of accepted confs of all processors.
520       not_done=.true.
521 C----------------------------------------------------------------------------
522 C Main loop.
523 c----------------------------------------------------------------------------
524       it=0
525       nout=0
526       LOOP1:do while (not_done)
527         it=it+1
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.
533         accepted=.false.
534         LOOP2:do while (.not. accepted)
535
536          LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
537           do i=1,nodenum-1
538            if(imtasks(i).eq.0) then
539             is=i
540             exit
541            endif
542           enddo
543 C Retrieve the angles of previously accepted conformation
544           do j=1,nvar
545             varia(j)=varold(j)
546           enddo
547           call var_to_geom(nvar,varia)
548 C Rebuild the chain.
549           call chainbuild
550 C Heat up the system, if necessary.
551           call heat(over)
552 C If temperature cannot be further increased, stop.
553           if (over) then 
554            finish=.true.
555           endif
556           MoveType=0
557           nbond=0
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
563            if (print_mc.gt.0)
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.'
571             finish=.true.
572             write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',
573      &         MoveType,' returned from PERTURB.'
574            endif
575            call chainbuild
576           else
577            MoveType=0
578            moves(0)=moves(0)+1
579            nstart_grow=iran_num(3,nres)
580            if (print_mc.gt.0)
581      &      write (iout,'(2a,i3)') 'Random-generated conformation',
582      &      ' - chain regrown from residue',nstart_grow
583            call gen_rand_conf(nstart_grow,*30)
584           endif
585           call geom_to_var(nvar,varia)
586           ngen=ngen+1
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
593            neneval=neneval+1 
594            accepted=.false.
595            noverlap=noverlap+1
596            if (noverlap.gt.maxoverlap) then
597             write (iout,*) 'Too many overlapping confs.',
598      &      ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
599             finish=.true.
600            endif
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,
605      &             CG_COMM, ierr)
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,
609      &             CG_COMM, ierr)
610 c           write (iout,*) 'it: tag=',tag
611            call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,
612      &             CG_COMM, ierr)
613 c           write (iout,*) 'eold: tag=',tag
614            call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, 
615      &             is, tag,
616      &             CG_COMM, ierr)
617 c           write (iout,*) 'varia: tag=',tag
618            call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, 
619      &             is, tag,
620      &             CG_COMM, ierr)
621 c           write (iout,*) 'varold: tag=',tag
622 #ifdef AIX
623            call flush_(iout)
624 #else
625            call flush(iout)
626 #endif
627            imtasks(is)=1
628            imtasks_n=imtasks_n+1
629 C End distribution
630           endif ! overlap
631          enddo LOOP3
632
633          flag = .false.
634          LOOP_RECV:do while(.not.flag)
635           do is=1, nodenum-1
636            call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
637            if(flag) then
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
659             imtasks(is)=0
660             imtasks_n=imtasks_n-1
661             exit 
662            endif
663           enddo
664          enddo LOOP_RECV
665
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,
673      &                      similar,EneLower)
674          if(iitt.ne.it.and..not.similar) then
675           call metropolis(nvar,varia,varold,etot,eold,accepted,
676      &                      similar,EneLower)
677           accepted=enelower
678          endif
679          if(etot.lt.eneglobal)eneglobal=etot
680 c         if(mod(it,100).eq.0)
681          write(iout,*)'CHUJOJEB ',neneval,eneglobal
682          if (accepted) then
683 C Write the accepted conformation.
684            nout=nout+1
685            if (refstr) then
686              call var_to_geom(nvar,varia)
687              call chainbuild
688              kkk=1
689              call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
690      &                    nsup,przes,obr,non_conv)
691              rms=dsqrt(rms)
692              call contact(.false.,ncont,icont,co)
693              frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
694              write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
695      &         'RMS deviation from the reference structure:',rms,
696      &         ' % of native contacts:',frac*100,' contact order:',co
697            endif ! refstr
698            if (print_mc.gt.0) 
699      &      write (iout,*) 'Writing new conformation',nout
700            if (print_stat) then
701              call var_to_geom(nvar,varia)
702 #if defined(AIX) || defined(PGI)
703              open (istat,file=statname,position='append')
704 #else
705              open (istat,file=statname,access='append')
706 #endif
707              if (refstr) then
708                write (istat,'(i5,16(1pe14.5))') nout,
709      &          (energia(print_order(i)),i=1,nprint_ene),
710      &          etot,rms,frac
711              else
712                write (istat,'(i5,16(1pe14.5))') nout,
713      &          (energia(print_order(i)),i=1,nprint_ene),etot
714              endif ! refstr
715              close(istat)
716            endif ! print_stat
717 C Print internal coordinates.
718            if (print_int) call briefout(nout,etot)
719            nacc=nacc+1
720            nacc_tot=nacc_tot+1
721            if (elowest.gt.etot) elowest=etot
722            moves_acc(MoveType)=moves_acc(MoveType)+1
723            if (MoveType.eq.1) then
724              nbond_acc(nbond)=nbond_acc(nbond)+1
725            endif
726 C Check against conformation repetitions.
727           irepet=conf_comp(varia,etot)
728           if (nrepm.gt.maxrepm) then
729            if (print_mc.gt.0) 
730      &      write (iout,'(a)') 'Too many conformation repetitions.'
731             finish=.true.
732            endif
733 C Store the accepted conf. and its energy.
734            eold=etot
735            do i=1,nvar
736             varold(i)=varia(i)
737            enddo
738            if (irepet.eq.0) call zapis(varia,etot)
739 C Lower the temperature, if necessary.
740            call cool
741           else
742            ntrial=ntrial+1
743          endif ! accepted
744    30    continue
745          if(finish.and.imtasks_n.eq.0)exit LOOP2
746         enddo LOOP2 ! accepted
747 C Check for time limit.
748         not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
749         if(.not.not_done .or. finish) then
750          if(imtasks_n.gt.0) then
751           not_done=.true.
752          else
753           not_done=.false.
754          endif
755          finish=.true.
756         endif
757       enddo LOOP1 ! not_done
758       runtime=tcpu()
759       if (print_mc.gt.0) then
760         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
761         call statprint(nacc,nfun,iretcode,etot,elowest)
762         write (iout,'(a)') 
763      & 'Statistics of multiple-bond motions. Total motions:' 
764         write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
765         write (iout,'(a)') 'Accepted motions:'
766         write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
767         if (it.ge.maxacc)
768      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
769       endif
770 #ifdef AIX
771       call flush_(iout)
772 #else
773       call flush(iout)
774 #endif
775       do is=1,nodenum-1
776         call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
777      &             CG_COMM, ierr)
778       enddo
779       return
780       end
781 c------------------------------------------------------------------------------
782       subroutine execute_slave(nodeinfo,iprint)
783       implicit real*8 (a-h,o-z)
784       include 'DIMENSIONS'
785       include 'mpif.h'
786       include 'COMMON.TIME1'
787       include 'COMMON.IOUNITS'
788 crc      include 'COMMON.DEFORM'
789 crc      include 'COMMON.DEFORM1'
790 crc      include 'COMMON.DEFORM2'
791       include 'COMMON.LOCAL'
792       include 'COMMON.VAR'
793       include 'COMMON.INFO'
794       include 'COMMON.MINIM'
795       character*10 nodeinfo 
796       double precision x(maxvar),x1(maxvar)
797       nodeinfo='chujwdupe'
798 c      print *,'Processor:',MyID,' Entering execute_slave'
799       tag=0
800 c      call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
801 c     &              CG_COMM, ierr)
802
803 1001  call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
804      &              CG_COMM, status, ierr)
805 c      write(iout,*)'12: tag=',tag
806       if(iprint.ge.2)print *, MyID,' recv order ',i_switch
807       if (i_switch.eq.12) then
808        i_grnum_d=0
809        i_ennum_d=0
810        i_hesnum_d=0
811        call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
812      &               CG_COMM, status, ierr)
813 c      write(iout,*)'12: tag=',tag
814        call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
815      &               CG_COMM, status, ierr)
816 c      write(iout,*)'ener: tag=',tag
817        call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
818      &               CG_COMM, status, ierr)
819 c      write(iout,*)'x: tag=',tag
820        call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
821      &               CG_COMM, status, ierr)
822 c      write(iout,*)'x1: tag=',tag
823 #ifdef AIX
824        call flush_(iout)
825 #else
826        call flush(iout)
827 #endif
828 c       print *,'calling minimize'
829        call minimize(energyx,x,iretcode,nfun)
830        if(iprint.gt.0)
831      &  write(iout,100)'minimized energy = ',energyx,
832      &    ' # funeval:',nfun,' iret ',iretcode
833         write(*,100)'minimized energy = ',energyx,
834      &    ' # funeval:',nfun,' iret ',iretcode
835  100   format(a20,f10.5,a12,i5,a6,i2)
836        if(iretcode.eq.10) then
837          do iminrep=2,3
838           if(iprint.gt.1)
839      &    write(iout,*)' ... not converged - trying again ',iminrep
840           call minimize(energyx,x,iretcode,nfun)
841           if(iprint.gt.1)
842      &    write(iout,*)'minimized energy = ',energyx,
843      &     ' # funeval:',nfun,' iret ',iretcode
844           if(iretcode.ne.10)go to 812
845          enddo
846          if(iretcode.eq.10) then
847           if(iprint.gt.1)
848      &    write(iout,*)' ... not converged again - giving up'
849           go to 812
850          endif
851        endif
852 812    continue
853 c       print *,'Sending results'
854        call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
855      &              CG_COMM, ierr)
856        call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
857      &              CG_COMM, ierr)
858        call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
859      &              CG_COMM, ierr)
860        call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
861      &              CG_COMM, ierr)
862        call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
863      &              CG_COMM, ierr)
864        call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
865      &              CG_COMM, ierr)
866        call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
867      &              CG_COMM, ierr)
868        call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
869      &              CG_COMM, ierr)
870 c       print *,'End sending'
871        go to 1001
872       endif
873
874       return
875       end
876 #endif
877 c------------------------------------------------------------------------------
878       subroutine statprint(it,nfun,iretcode,etot,elowest)
879       implicit real*8 (a-h,o-z)
880       include 'DIMENSIONS'
881       include 'COMMON.IOUNITS'
882       include 'COMMON.CONTROL'
883       include 'COMMON.MCM'
884       if (minim) then
885         write (iout,
886      &  '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
887      &  'Finished iteration #',it,' energy is',etot,
888      &  ' lowest energy:',elowest,
889      &  'SUMSL return code:',iretcode,
890      &  ' # of energy evaluations:',neneval,
891      &  '# of temperature jumps:',ntherm,
892      &  ' # of minima repetitions:',nrepm
893       else
894         write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
895      &  'Finished iteration #',it,' energy is',etot,
896      &  ' lowest energy:',elowest
897       endif
898       write (iout,'(/4a)')
899      & 'Kind of move   ','           total','       accepted',
900      & '  fraction'
901       write (iout,'(58(1h-))')
902       do i=-1,MaxMoveType
903         if (moves(i).eq.0) then
904           fr_mov_i=0.0d0
905         else
906           fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
907         endif
908         write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
909      &         fr_mov_i
910       enddo
911       write (iout,'(a,2i15,f10.5)') 'total           ',nmove,nacc_tot,
912      &         dfloat(nacc_tot)/dfloat(nmove)
913       write (iout,'(58(1h-))')
914       write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
915       return
916       end
917 c------------------------------------------------------------------------------
918       subroutine heat(over)
919       implicit real*8 (a-h,o-z)
920       include 'DIMENSIONS'
921       include 'COMMON.MCM'
922       include 'COMMON.IOUNITS'
923       logical over
924 C Check if there`s a need to increase temperature.
925       if (ntrial.gt.maxtrial) then
926         if (NstepH.gt.0) then
927           if (dabs(Tcur-TMax).lt.1.0D-7) then
928             if (print_mc.gt.0)
929      &      write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') 
930      &      'Upper limit of temperature reached. Terminating.'
931             over=.true.
932             Tcur=Tmin
933           else
934             Tcur=Tcur*TstepH
935             if (Tcur.gt.Tmax) Tcur=Tmax
936             betbol=1.0D0/(Rbol*Tcur)
937             if (print_mc.gt.0)
938      &      write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
939      &      'System heated up to ',Tcur,' K; BetBol:',betbol
940             ntherm=ntherm+1
941             ntrial=0
942             over=.false.
943           endif
944         else
945          if (print_mc.gt.0)
946      &   write (iout,'(a)') 
947      & 'Maximum number of trials in a single MCM iteration exceeded.'
948          over=.true.
949          Tcur=Tmin
950         endif
951       else
952         over=.false.
953       endif
954       return
955       end
956 c------------------------------------------------------------------------------
957       subroutine cool
958       implicit real*8 (a-h,o-z)
959       include 'DIMENSIONS'
960       include 'COMMON.MCM'
961       include 'COMMON.IOUNITS'
962       if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
963         Tcur=Tcur/TstepC
964         if (Tcur.lt.Tmin) Tcur=Tmin
965         betbol=1.0D0/(Rbol*Tcur)
966         if (print_mc.gt.0)
967      &  write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
968      &  'System cooled down up to ',Tcur,' K; BetBol:',betbol
969       endif
970       return
971       end
972 C---------------------------------------------------------------------------
973       subroutine zapis(varia,etot)
974       implicit real*8 (a-h,o-z)
975       include 'DIMENSIONS'
976 #ifdef MP
977       include 'mpif.h'
978       include 'COMMON.INFO'
979 #endif
980       include 'COMMON.GEO'
981       include 'COMMON.VAR'
982       include 'COMMON.MCM'
983       include 'COMMON.IOUNITS'
984       integer itemp(maxsave)
985       double precision varia(maxvar)
986       logical lprint
987       lprint=.false.
988       if (lprint) then
989       write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
990      &  ' MaxSave=',MaxSave
991       write (iout,'(a)') 'Current energy and conformation:'
992       write (iout,'(1pe14.5)') etot
993       write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
994       endif
995 C Shift the contents of the esave and varsave arrays if filled up.
996       call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
997      &               etot,varia,esave,varsave)
998       if (lprint) then
999       write (iout,'(a)') 'Energies and the VarSave array.'
1000       do i=1,nsave
1001         write (iout,'(i5,1pe14.5)') i,esave(i)
1002         write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
1003       enddo
1004       endif
1005       return
1006       end 
1007 C---------------------------------------------------------------------------
1008       subroutine perturb(error,lprint,MoveType,nbond,max_phi)
1009       implicit real*8 (a-h,o-z)
1010       include 'DIMENSIONS'
1011       parameter (MMaxSideMove=100)
1012       include 'COMMON.MCM'
1013       include 'COMMON.CHAIN'
1014       include 'COMMON.INTERACT'
1015       include 'COMMON.VAR'
1016       include 'COMMON.GEO'
1017       include 'COMMON.NAMES'
1018       include 'COMMON.IOUNITS'
1019 crc      include 'COMMON.DEFORM1'
1020       logical error,lprint,fail
1021       integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
1022       double precision max_phi
1023       double precision psi,gen_psi
1024       external iran_num
1025       integer iran_num
1026       integer ifour 
1027       data ifour /4/
1028       error=.false.
1029       lprint=.false.
1030 C Perturb the conformation according to a randomly selected move.
1031       call SelectMove(MoveType)
1032 c      write (iout,*) 'MoveType=',MoveType
1033       itrial=0
1034       goto (100,200,300,400,500) MoveType
1035 C------------------------------------------------------------------------------
1036 C Backbone N-bond move.
1037 C Select the number of bonds (length of the segment to perturb).
1038   100 continue
1039       if (itrial.gt.1000) then
1040         write (iout,'(a)') 'Too many attempts at multiple-bond move.'
1041         error=.true.
1042         return
1043       endif
1044       bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
1045 c     print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
1046 c    & ' Bond_prob=',Bond_Prob
1047       do i=1,nbm-1
1048 c       print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
1049         if (bond_prob.ge.sumpro_bond(i) .and. 
1050      &               bond_prob.le.sumpro_bond(i+1)) then
1051           nbond=i+1
1052           goto 10
1053         endif
1054       enddo
1055       write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
1056      &                    ' to move out of range.'
1057       error=.true.
1058       return
1059    10 continue
1060       if (nwindow.gt.0) then
1061 C Select the first residue to perturb
1062         iwindow=iran_num(1,nwindow)
1063         print *,'iwindow=',iwindow
1064         iiwin=1
1065         do while (winlen(iwindow).lt.nbond)
1066           iwindow=iran_num(1,nwindow)
1067           iiwin=iiwin+1
1068           if (iiwin.gt.1000) then
1069              write (iout,'(a)') 'Cannot select moveable residues.'
1070              error=.true.
1071              return
1072           endif
1073         enddo 
1074         nstart=iran_num(winstart(iwindow),winend(iwindow))
1075       else
1076         nstart = iran_num(koniecl+2,nres-nbond-koniecl)  
1077 cd      print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
1078 cd   &        ' nstart=',nstart
1079       endif
1080       psi = gen_psi()
1081       if (psi.eq.0.0) then
1082         error=.true.
1083         return
1084       endif
1085       if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
1086      & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
1087 cd    print *,'nstart=',nstart
1088       call bond_move(nbond,nstart,psi,.false.,error)
1089       if (error) then 
1090         write (iout,'(2a)') 
1091      & 'Could not define reference system in bond_move, ',
1092      & 'choosing ahother segment.'
1093         itrial=itrial+1
1094         goto 100
1095       endif
1096       nbond_move(nbond)=nbond_move(nbond)+1
1097       moves(1)=moves(1)+1
1098       nmove=nmove+1
1099       return
1100 C------------------------------------------------------------------------------
1101 C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
1102 C the chain.
1103   200 continue
1104       lprint=.true.
1105 c     end_select=iran_num(1,2*koniecl)
1106 c     if (end_select.gt.koniecl) then
1107 c       end_select=nphi-(end_select-koniecl)
1108 c     else 
1109 c       end_select=koniecl+3
1110 c     endif
1111 c     if (nwindow.gt.0) then
1112 c       iwin=iran_num(1,nwindow)
1113 c       i1=max0(4,winstart(iwin))
1114 c       i2=min0(winend(imin)+2,nres)
1115 c       end_select=iran_num(i1,i2)
1116 c     else
1117 c      iselect = iran_num(1,nmov_var)
1118 c      jj = 0
1119 c      do i=1,nphi
1120 c        if (isearch_tab(i).eq.1) jj = jj+1
1121 c        if (jj.eq.iselect) then
1122 c          end_select=i+3
1123 c          exit
1124 c        endif
1125 c      enddo    
1126 c     endif
1127       end_select = iran_num(4,nres)
1128       psi=max_phi*gen_psi()
1129       if (psi.eq.0.0D0) then
1130         error=.true.
1131         return
1132       endif
1133       phi(end_select)=pinorm(phi(end_select)+psi)
1134       if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') 
1135      & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
1136      & phi(end_select)*rad2deg   
1137 c     if (end_select.gt.3) 
1138 c    &   theta(end_select-1)=gen_theta(itype(end_select-2),
1139 c    &                              phi(end_select-1),phi(end_select))
1140 c     if (end_select.lt.nres) 
1141 c    &    theta(end_select)=gen_theta(itype(end_select-1),
1142 c    &                              phi(end_select),phi(end_select+1))
1143 cd    print *,'nres=',nres,' end_select=',end_select
1144 cd    print *,'theta',end_select-1,theta(end_select-1)
1145 cd    print *,'theta',end_select,theta(end_select)
1146       moves(2)=moves(2)+1
1147       nmove=nmove+1
1148       lprint=.false.
1149       return
1150 C------------------------------------------------------------------------------
1151 C Side chain move.
1152 C Select the number of SCs to perturb.
1153   300 isctry=0 
1154   301 nside_move=iran_num(1,MaxSideMove) 
1155 c     print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
1156 C Select the indices.
1157       do i=1,nside_move
1158         icount=0
1159   111   inds=iran_num(nnt,nct) 
1160         icount=icount+1
1161         if (icount.gt.1000) then
1162           write (iout,'(a)')'Error - cannot select side chains to move.'
1163           error=.true.
1164           return
1165         endif
1166         if (itype(inds).eq.10) goto 111
1167         do j=1,i-1
1168           if (inds.eq.ind_side(j)) goto 111
1169         enddo
1170         do j=1,i-1
1171           if (inds.lt.ind_side(j)) then
1172             indx=j
1173             goto 112
1174           endif
1175         enddo
1176         indx=i
1177   112   do j=i,indx+1,-1
1178           ind_side(j)=ind_side(j-1)
1179         enddo 
1180   113   ind_side(indx)=inds
1181       enddo
1182 C Carry out perturbation.
1183       do i=1,nside_move
1184         ii=ind_side(i)
1185         iti=itype(ii)
1186         call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
1187         if (fail) then
1188           isctry=isctry+1
1189           if (isctry.gt.1000) then
1190             write (iout,'(a)') 'Too many errors in SC generation.'
1191             error=.true.
1192             return
1193           endif
1194           goto 301 
1195         endif
1196         if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') 
1197      &   'Side chain ',restyp(iti),ii,' moved to ',
1198      &   alph(ii)*rad2deg,omeg(ii)*rad2deg
1199       enddo
1200       moves(3)=moves(3)+1
1201       nmove=nmove+1
1202       return
1203 C------------------------------------------------------------------------------
1204 C THETA move
1205   400 end_select=iran_num(3,nres)
1206       theta_new=gen_theta(itype(end_select),phi(end_select),
1207      &                    phi(end_select+1))
1208       if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') 
1209      & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
1210      & ' to ',theta_new*rad2deg
1211       theta(end_select)=theta_new  
1212       moves(4)=moves(4)+1
1213       nmove=nmove+1 
1214       return
1215 C------------------------------------------------------------------------------
1216 C Error returned from SelectMove.
1217   500 error=.true.
1218       return
1219       end
1220 C------------------------------------------------------------------------------
1221       subroutine SelectMove(MoveType)
1222       implicit real*8 (a-h,o-z)
1223       include 'DIMENSIONS'
1224       include 'COMMON.MCM'
1225       include 'COMMON.IOUNITS'
1226       what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
1227       do i=1,MaxMoveType
1228         if (what_move.ge.sumpro_type(i-1).and.
1229      &            what_move.lt.sumpro_type(i)) then
1230           MoveType=i
1231           return
1232         endif
1233       enddo
1234       write (iout,'(a)') 
1235      & 'Fatal error in SelectMoveType: cannot select move.'
1236       MoveType=MaxMoveType+1
1237       return
1238       end
1239 c----------------------------------------------------------------------------
1240       double precision function gen_psi()
1241       implicit none
1242       integer i
1243       double precision x,ran_number
1244       include 'COMMON.IOUNITS'
1245       include 'COMMON.GEO'
1246       x=0.0D0
1247       do i=1,100
1248         x=ran_number(-pi,pi)
1249         if (dabs(x).gt.angmin) then
1250           gen_psi=x
1251           return
1252         endif
1253       enddo
1254       write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
1255       gen_psi=0.0D0
1256       return
1257       end
1258 c----------------------------------------------------------------------------
1259       subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
1260      &                      enelower)
1261       implicit real*8 (a-h,o-z)
1262       include 'DIMENSIONS'
1263       include 'COMMON.MCM'
1264       include 'COMMON.IOUNITS'
1265       include 'COMMON.VAR'
1266       include 'COMMON.GEO'
1267 crc      include 'COMMON.DEFORM'
1268       double precision ecur,eold,xx,ran_number,bol
1269       double precision xcur(n),xold(n)
1270       double precision ecut1 ,ecut2 ,tola
1271       logical accepted,similar,not_done,enelower
1272       logical lprn
1273       data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
1274 !      ecut1=-5*enedif
1275 !      ecut2=50*enedif
1276 !      tola=5.0d0
1277 C Set lprn=.true. for debugging.
1278       lprn=.false.
1279       if (lprn) 
1280      &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
1281       similar=.false.
1282       enelower=.false.
1283       accepted=.false.
1284 C Check if the conformation is similar.
1285       difene=ecur-eold
1286       reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
1287       if (lprn) then
1288         write (iout,*) 'Metropolis'
1289         write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
1290       endif
1291 C If energy went down remarkably, we accept the new conformation 
1292 C unconditionally.
1293 cjp      if (reldife.lt.ecut1) then
1294       if (difene.lt.ecut1) then
1295         accepted=.true.
1296         EneLower=.true.
1297         if (lprn) write (iout,'(a)') 
1298      &   'Conformation accepted, because energy has lowered remarkably.'
1299 !      elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola) 
1300 cjp      elseif (reldife.lt.ecut2) 
1301       elseif (difene.lt.ecut2) 
1302      & then
1303 C Reject the conf. if energy has changed insignificantly and there is not 
1304 C much change in conformation.
1305         if (lprn) 
1306      &   write (iout,'(2a)') 'Conformation rejected, because it is',
1307      &      ' similar to the preceding one.'
1308         accepted=.false.
1309         similar=.true.
1310       else 
1311 C Else carry out Metropolis test.
1312         EneLower=.false.
1313         xx=ran_number(0.0D0,1.0D0) 
1314         xxh=betbol*difene
1315         if (lprn)
1316      &    write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
1317         if (xxh.gt.50.0D0) then
1318           bol=0.0D0
1319         else
1320           bol=exp(-xxh)
1321         endif
1322         if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
1323         if (bol.gt.xx) then
1324           accepted=.true. 
1325           if (lprn) write (iout,'(a)') 
1326      &    'Conformation accepted, because it passed Metropolis test.'
1327         else
1328           accepted=.false.
1329           if (lprn) write (iout,'(a)') 
1330      & 'Conformation rejected, because it did not pass Metropolis test.'
1331         endif
1332       endif
1333 #ifdef AIX
1334       call flush_(iout)
1335 #else
1336       call flush(iout)
1337 #endif
1338       return
1339       end 
1340 c------------------------------------------------------------------------------
1341       integer function conf_comp(x,ene)
1342       implicit real*8 (a-h,o-z)
1343       include 'DIMENSIONS'
1344       include 'COMMON.MCM'
1345       include 'COMMON.VAR'
1346       include 'COMMON.IOUNITS'
1347       include 'COMMON.GEO' 
1348       double precision etol , angtol 
1349       double precision x(maxvar)
1350       double precision dif_ang,difa
1351       data etol /0.1D0/, angtol /20.0D0/
1352       do ii=nsave,1,-1
1353 c       write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
1354         if (dabs(ene-esave(ii)).lt.etol) then
1355           difa=dif_ang(nphi,x,varsave(1,ii))
1356 c         do i=1,nphi
1357 c           write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1358 c    &          rad2deg*varsave(i,ii)
1359 c         enddo
1360 c         write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
1361           if (difa.le.angtol) then
1362             if (print_mc.gt.0) then
1363             write (iout,'(a,i5,2(a,1pe15.4))') 
1364      &      'Current conformation matches #',ii,
1365      &      ' in the store array ene=',ene,' esave=',esave(ii)
1366 c           write (*,'(a,i5,a)') 'Current conformation matches #',ii,
1367 c    &      ' in the store array.'
1368             endif ! print_mc.gt.0
1369             if (print_mc.gt.1) then
1370             do i=1,nphi
1371               write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1372      &            rad2deg*varsave(i,ii)
1373             enddo
1374             endif ! print_mc.gt.1
1375             nrepm=nrepm+1
1376             conf_comp=ii
1377             return
1378           endif
1379         endif
1380       enddo 
1381       conf_comp=0
1382       return
1383       end 
1384 C----------------------------------------------------------------------------
1385       double precision function dif_ang(n,x,y)
1386       implicit none
1387       integer i,n
1388       double precision x(n),y(n)
1389       double precision w,wa,dif,difa
1390       double precision pinorm 
1391       include 'COMMON.GEO'
1392       wa=0.0D0
1393       difa=0.0D0
1394       do i=1,n
1395         dif=dabs(pinorm(y(i)-x(i)))
1396         if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
1397         w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
1398         wa=wa+w
1399         difa=difa+dif*dif*w
1400       enddo 
1401       dif_ang=rad2deg*dsqrt(difa/wa)
1402       return
1403       end
1404 c--------------------------------------------------------------------------
1405       subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
1406      &                     ecur,xcur,ecache,xcache)
1407       implicit none 
1408       include 'COMMON.GEO'
1409       include 'COMMON.IOUNITS'
1410       integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
1411       integer i,ii,j
1412       double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2) 
1413 cd    write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
1414 cd    write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
1415 cd    write (iout,*) 'Old CACHE array:'
1416 cd    do i=1,ncache
1417 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1418 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1419 cd    enddo
1420
1421       i=ncache
1422       do while (i.gt.0 .and. ecur.lt.ecache(i))
1423         i=i-1
1424       enddo
1425       i=i+1
1426 cd    write (iout,*) 'i=',i,' ncache=',ncache
1427       if (ncache.eq.n2) then
1428         write (iout,*) 'Cache dimension exceeded',ncache,n2
1429         write (iout,*) 'Highest-energy conformation will be removed.'
1430         ncache=ncache-1
1431       endif
1432       do ii=ncache,i,-1
1433         ecache(ii+1)=ecache(ii)
1434         CachSrc(ii+1)=CachSrc(ii)
1435         do j=1,nvar
1436           xcache(j,ii+1)=xcache(j,ii)
1437         enddo
1438       enddo
1439       ecache(i)=ecur
1440       CachSrc(i)=SourceID
1441       do j=1,nvar
1442         xcache(j,i)=xcur(j)
1443       enddo
1444       ncache=ncache+1
1445 cd    write (iout,*) 'New CACHE array:'
1446 cd    do i=1,ncache
1447 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1448 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1449 cd    enddo
1450       return
1451       end
1452 c--------------------------------------------------------------------------
1453       subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
1454      &                         xcache)
1455       implicit none 
1456       include 'COMMON.GEO'
1457       include 'COMMON.IOUNITS'
1458       integer n1,n2,ncache,nvar,CachSrc(n2)
1459       integer i,ii,j
1460       double precision ecache(n2),xcache(n1,n2) 
1461
1462 cd    write (iout,*) 'Enter RM_FROM_CACHE'
1463 cd    write (iout,*) 'Old CACHE array:'
1464 cd    do ii=1,ncache
1465 cd    write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
1466 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
1467 cd    enddo
1468
1469       do ii=i+1,ncache
1470         ecache(ii-1)=ecache(ii)
1471         CachSrc(ii-1)=CachSrc(ii)
1472         do j=1,nvar
1473           xcache(j,ii-1)=xcache(j,ii)
1474         enddo
1475       enddo
1476       ncache=ncache-1
1477 cd    write (iout,*) 'New CACHE array:'
1478 cd    do i=1,ncache
1479 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1480 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1481 cd    enddo
1482       return
1483       end