Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / unres / src_MD-M / unres.F
1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2 C                                                                              C
3 C                                U N R E S                                     C
4 C                                                                              C
5 C Program to carry out conformational search of proteins in an united-residue  C
6 C approximation.                                                               C
7 C                                                                              C
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9       implicit real*8 (a-h,o-z)
10       include 'DIMENSIONS'
11
12
13 #ifdef MPI
14       include 'mpif.h'
15       include 'COMMON.SETUP'
16 #endif
17       include 'COMMON.TIME1'
18       include 'COMMON.INTERACT'
19       include 'COMMON.NAMES'
20       include 'COMMON.GEO'
21       include 'COMMON.HEADER'
22       include 'COMMON.CONTROL'
23       include 'COMMON.CONTACTS'
24       include 'COMMON.CHAIN'
25       include 'COMMON.VAR'
26       include 'COMMON.IOUNITS'
27       include 'COMMON.FFIELD'
28       include 'COMMON.REMD'
29       include 'COMMON.MD'
30       include 'COMMON.SBRIDGE'
31       double precision hrtime,mintime,sectime
32       character*64 text_mode_calc(-2:14) /'test',
33      & 'SC rotamer distribution',
34      & 'Energy evaluation or minimization',
35      & 'Regularization of PDB structure',
36      & 'Threading of a sequence on PDB structures',
37      & 'Monte Carlo (with minimization) ',
38      & 'Energy minimization of multiple conformations',
39      & 'Checking energy gradient',
40      & 'Entropic sampling Monte Carlo (with minimization)',
41      & 'Energy map',
42      & 'CSA calculations',
43      & 'Not used 9',
44      & 'Not used 10',
45      & 'Soft regularization of PDB structure',
46      & 'Mesoscopic molecular dynamics (MD) ',
47      & 'Not used 13',
48      & 'Replica exchange molecular dynamics (REMD)'/
49       external ilen
50
51 c      call memmon_print_usage()
52
53       call init_task
54       if (me.eq.king)
55      & write(iout,*)'### LAST MODIFIED  4/25/08 7:29PM by adam'  
56       if (me.eq.king) call cinfo
57 C Read force field parameters and job setup data
58       call readrtns
59       call flush(iout)
60 C
61       if (me.eq.king .or. .not. out1file) then
62        write (iout,'(2a/)') 
63      & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
64      & ' calculation.' 
65        if (minim) write (iout,'(a)') 
66      &  'Conformations will be energy-minimized.'
67        write (iout,'(80(1h*)/)') 
68       endif
69       call flush(iout)
70 C
71       if (modecalc.eq.-2) then
72         call test
73         stop
74       else if (modecalc.eq.-1) then
75         write(iout,*) "call check_sc_map next"
76         call check_bond
77         stop
78       endif
79 #ifdef MPI
80       if (fg_rank.gt.0) then
81 C Fine-grain slaves just do energy and gradient components.
82         call ergastulum ! slave workhouse in Latin
83       else
84 #endif
85       if (modecalc.eq.0) then
86         call exec_eeval_or_minim
87       else if (modecalc.eq.1) then
88         call exec_regularize
89       else if (modecalc.eq.2) then
90         call exec_thread
91       else if (modecalc.eq.3 .or. modecalc .eq.6) then
92         call exec_MC
93       else if (modecalc.eq.4) then
94         call exec_mult_eeval_or_minim
95       else if (modecalc.eq.5) then
96         call exec_checkgrad
97       else if (ModeCalc.eq.7) then
98         call exec_map
99       else if (ModeCalc.eq.8) then
100         call exec_CSA
101       else if (modecalc.eq.11) then
102         call exec_softreg
103       else if (modecalc.eq.12) then
104         call exec_MD
105       else if (modecalc.eq.14) then
106         call exec_MREMD
107       else
108         write (iout,'(a)') 'This calculation type is not supported',
109      &   ModeCalc
110       endif
111 #ifdef MPI
112       endif
113 C Finish task.
114       if (fg_rank.eq.0) call finish_task
115 c      call memmon_print_usage()
116 #ifdef TIMING
117        call print_detailed_timing
118 #endif
119       call MPI_Finalize(ierr)
120       stop 'Bye Bye...'
121 #else
122       call dajczas(tcpu(),hrtime,mintime,sectime)
123       stop '********** Program terminated normally.'
124 #endif
125       end
126 c--------------------------------------------------------------------------
127       subroutine exec_MD
128       include 'DIMENSIONS'
129 #ifdef MPI
130       include "mpif.h"
131 #endif
132       include 'COMMON.SETUP'
133       include 'COMMON.CONTROL'
134       include 'COMMON.IOUNITS'
135       if (me.eq.king .or. .not. out1file)
136      &   write (iout,*) "Calling chainbuild"
137       call chainbuild
138       call MD
139       return
140       end
141 c---------------------------------------------------------------------------
142       subroutine exec_MREMD
143       include 'DIMENSIONS'
144 #ifdef MPI
145       include "mpif.h"
146 #endif
147       include 'COMMON.SETUP'
148       include 'COMMON.CONTROL'
149       include 'COMMON.IOUNITS'
150       include 'COMMON.REMD'
151       if (me.eq.king .or. .not. out1file)
152      &   write (iout,*) "Calling chainbuild"
153       call chainbuild
154       if (me.eq.king .or. .not. out1file)
155      &   write (iout,*) "Calling REMD"
156       if (remd_mlist) then 
157         call MREMD
158       else
159         do i=1,nrep
160           remd_m(i)=1
161         enddo
162         call MREMD
163       endif
164       return
165       end
166 c---------------------------------------------------------------------------
167       subroutine exec_eeval_or_minim
168       implicit real*8 (a-h,o-z)
169       include 'DIMENSIONS'
170 #ifdef MPI
171       include 'mpif.h'
172 #endif
173       include 'COMMON.SETUP'
174       include 'COMMON.TIME1'
175       include 'COMMON.INTERACT'
176       include 'COMMON.NAMES'
177       include 'COMMON.GEO'
178       include 'COMMON.HEADER'
179       include 'COMMON.CONTROL'
180       include 'COMMON.CONTACTS'
181       include 'COMMON.CHAIN'
182       include 'COMMON.VAR'
183       include 'COMMON.IOUNITS'
184       include 'COMMON.FFIELD'
185       include 'COMMON.REMD'
186       include 'COMMON.MD'
187       include 'COMMON.SBRIDGE'
188       common /srutu/ icall
189       double precision energy(0:n_ene)
190       double precision energy_long(0:n_ene),energy_short(0:n_ene)
191       double precision varia(maxvar)
192       if (indpdb.eq.0) call chainbuild
193       time00=MPI_Wtime()
194       call chainbuild_cart
195       if (split_ene) then
196        print *,"Processor",myrank," after chainbuild"
197        icall=1
198        call etotal_long(energy_long(0))
199        write (iout,*) "Printing long range energy"
200        call enerprint(energy_long(0))
201        call etotal_short(energy_short(0))
202        write (iout,*) "Printing short range energy"
203        call enerprint(energy_short(0))
204        do i=0,n_ene
205          energy(i)=energy_long(i)+energy_short(i)
206          write (iout,*) i,energy_long(i),energy_short(i),energy(i)
207        enddo
208        write (iout,*) "Printing long+short range energy"
209        call enerprint(energy(0))
210       endif
211       call etotal(energy(0))
212       time_ene=MPI_Wtime()-time00
213       write (iout,*) "Time for energy evaluation",time_ene
214       print *,"after etotal"
215       etota = energy(0)
216       etot =etota
217       call enerprint(energy(0))
218       call hairpin(.true.,nharp,iharp)
219       call secondary2(.true.)
220       if (minim) then
221 crc overlap test
222         if (overlapsc) then 
223           print *, 'Calling OVERLAP_SC'
224           call overlap_sc(fail)
225         endif 
226
227         if (searchsc) then 
228           call sc_move(2,nres-1,10,1d10,nft_sc,etot)
229           print *,'SC_move',nft_sc,etot
230           write(iout,*) 'SC_move',nft_sc,etot
231         endif 
232
233         if (dccart) then
234           print *, 'Calling MINIM_DC'
235           time1=MPI_WTIME()
236           call minim_dc(etot,iretcode,nfun)
237         else
238           if (indpdb.ne.0) then 
239             call bond_regular
240             call chainbuild
241           endif
242           call geom_to_var(nvar,varia)
243           print *,'Calling MINIMIZE.'
244           time1=MPI_WTIME()
245           call minimize(etot,varia,iretcode,nfun)
246         endif
247         print *,'SUMSL return code is',iretcode,' eval ',nfun
248         evals=nfun/(MPI_WTIME()-time1)
249         print *,'# eval/s',evals
250         print *,'refstr=',refstr
251         call hairpin(.true.,nharp,iharp)
252         call secondary2(.true.)
253         call etotal(energy(0))
254         etot = energy(0)
255         call enerprint(energy(0))
256
257         call intout
258         call briefout(0,etot)
259         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
260           write (iout,'(a,i3)') 'SUMSL return code:',iretcode
261           write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
262           write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
263       else
264         print *,'refstr=',refstr
265         if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
266         call briefout(0,etot)
267       endif
268       if (outpdb) call pdbout(etot,titel(:32),ipdb)
269       if (outmol2) call mol2out(etot,titel(:32))
270       return
271       end
272 c---------------------------------------------------------------------------
273       subroutine exec_regularize
274       implicit real*8 (a-h,o-z)
275       include 'DIMENSIONS'
276 #ifdef MPI
277       include 'mpif.h'
278 #endif
279       include 'COMMON.SETUP'
280       include 'COMMON.TIME1'
281       include 'COMMON.INTERACT'
282       include 'COMMON.NAMES'
283       include 'COMMON.GEO'
284       include 'COMMON.HEADER'
285       include 'COMMON.CONTROL'
286       include 'COMMON.CONTACTS'
287       include 'COMMON.CHAIN'
288       include 'COMMON.VAR'
289       include 'COMMON.IOUNITS'
290       include 'COMMON.FFIELD'
291       include 'COMMON.REMD'
292       include 'COMMON.MD'
293       include 'COMMON.SBRIDGE'
294       double precision energy(0:n_ene)
295
296       call gen_dist_constr
297       call sc_conf
298       call intout
299       call regularize(nct-nnt+1,etot,rms,cref(1,nnt,1),iretcode)
300       call etotal(energy(0))
301       energy(0)=energy(0)-energy(14)
302       etot=energy(0)
303       call enerprint(energy(0))
304       call intout
305       call briefout(0,etot)
306       if (outpdb) call pdbout(etot,titel(:32),ipdb)
307       if (outmol2) call mol2out(etot,titel(:32))
308       if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
309       write (iout,'(a,i3)') 'SUMSL return code:',iretcode
310       return
311       end
312 c---------------------------------------------------------------------------
313       subroutine exec_thread
314       include 'DIMENSIONS'
315 #ifdef MP
316       include "mpif.h"
317 #endif
318       include "COMMON.SETUP"
319       call thread_seq
320       return
321       end
322 c---------------------------------------------------------------------------
323       subroutine exec_MC
324       implicit real*8 (a-h,o-z)
325       include 'DIMENSIONS'
326       character*10 nodeinfo
327       double precision varia(maxvar)
328 #ifdef MPI
329       include "mpif.h"
330 #endif
331       include "COMMON.SETUP"
332       include 'COMMON.CONTROL'
333       call mcm_setup
334       if (minim) then
335 #ifdef MPI
336         if (modecalc.eq.3) then
337           call do_mcm(ipar)
338         else
339           call entmcm
340         endif
341 #else
342         if (modecalc.eq.3) then
343           call do_mcm(ipar)
344         else
345           call entmcm
346         endif
347 #endif
348       else
349         call monte_carlo
350       endif
351       return
352       end
353 c---------------------------------------------------------------------------
354       subroutine exec_mult_eeval_or_minim
355       implicit real*8 (a-h,o-z)
356       include 'DIMENSIONS'
357 #ifdef MPI
358       include 'mpif.h'
359       dimension muster(mpi_status_size)
360 #endif
361       include 'COMMON.SETUP'
362       include 'COMMON.TIME1'
363       include 'COMMON.INTERACT'
364       include 'COMMON.NAMES'
365       include 'COMMON.GEO'
366       include 'COMMON.HEADER'
367       include 'COMMON.CONTROL'
368       include 'COMMON.CONTACTS'
369       include 'COMMON.CHAIN'
370       include 'COMMON.VAR'
371       include 'COMMON.IOUNITS'
372       include 'COMMON.FFIELD'
373       include 'COMMON.REMD'
374       include 'COMMON.MD'
375       include 'COMMON.SBRIDGE'
376       double precision varia(maxvar)
377       dimension ind(6)
378       double precision energy(0:max_ene)
379       logical eof
380       eof=.false.
381 #ifdef MPI
382       if(me.ne.king) then
383         call minim_mcmf
384         return
385       endif
386
387       close (intin)
388       open(intin,file=intinname,status='old')
389       write (istat,'(a5,20a12)')"#    ",
390      &  (wname(print_order(i)),i=1,nprint_ene)
391       if (refstr) then
392         write (istat,'(a5,20a12)')"#    ",
393      &   (ename(print_order(i)),i=1,nprint_ene),
394      &   "ETOT total","RMSD","nat.contact","nnt.contact"        
395       else
396         write (istat,'(a5,20a12)')"#    ",
397      &    (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
398       endif
399
400       if (.not.minim) then
401         do while (.not. eof)
402           if (read_cart) then
403             read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
404             call read_x(intin,*11)
405 #ifdef MPI
406 c Broadcast the order to compute internal coordinates to the slaves.
407             if (nfgtasks.gt.1)
408      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
409 #endif
410             call int_from_cart1(.false.)
411           else
412             read (intin,'(i5)',end=1100,err=1100) iconf
413             call read_angles(intin,*11)
414             call geom_to_var(nvar,varia)
415             call chainbuild
416           endif
417           write (iout,'(a,i7)') 'Conformation #',iconf
418           call etotal(energy(0))
419           call briefout(iconf,energy(0))
420           call enerprint(energy(0))
421           etot=energy(0)
422           if (refstr) then 
423             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
424             write (istat,'(i5,20(f12.3))') iconf,
425      &      (energy(print_order(i)),i=1,nprint_ene),etot,
426      &       rms,frac,frac_nn,co
427 cjlee end
428           else
429             write (istat,'(i5,16(f12.3))') iconf,
430      &     (energy(print_order(i)),i=1,nprint_ene),etot
431           endif
432         enddo
433 1100    continue
434         goto 1101
435       endif
436
437       mm=0
438       imm=0
439       nft=0
440       ene0=0.0d0
441       n=0
442       iconf=0
443 c      do n=1,nzsc
444       do while (.not. eof)
445         mm=mm+1
446         if (mm.lt.nodes) then
447           if (read_cart) then
448             read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
449             call read_x(intin,*11)
450 #ifdef MPI
451 c Broadcast the order to compute internal coordinates to the slaves.
452             if (nfgtasks.gt.1) 
453      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
454 #endif
455             call int_from_cart1(.false.)
456           else
457             read (intin,'(i5)',end=11,err=11) iconf
458             call read_angles(intin,*11)
459             call geom_to_var(nvar,varia)
460             call chainbuild
461           endif
462           write (iout,'(a,i7)') 'Conformation #',iconf
463           n=n+1
464          imm=imm+1
465          ind(1)=1
466          ind(2)=n
467          ind(3)=0
468          ind(4)=0
469          ind(5)=0
470          ind(6)=0
471          ene0=0.0d0
472          call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
473      *                  ierr)
474          call mpi_send(varia,nvar,mpi_double_precision,mm,
475      *                  idreal,CG_COMM,ierr)
476          call mpi_send(ene0,1,mpi_double_precision,mm,
477      *                  idreal,CG_COMM,ierr)
478 c         print *,'task ',n,' sent to worker ',mm,nvar
479         else
480          call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
481      *                 CG_COMM,muster,ierr)
482          man=muster(mpi_source)
483 c         print *,'receiving result from worker ',man,' (',iii1,iii,')'
484          call mpi_recv(varia,nvar,mpi_double_precision, 
485      *               man,idreal,CG_COMM,muster,ierr)
486          call mpi_recv(ene,1,
487      *               mpi_double_precision,man,idreal,
488      *               CG_COMM,muster,ierr)
489          call mpi_recv(ene0,1,
490      *               mpi_double_precision,man,idreal,
491      *               CG_COMM,muster,ierr)
492 c         print *,'result received from worker ',man,' sending now'
493
494           call var_to_geom(nvar,varia)
495           call chainbuild
496           call etotal(energy(0))
497           iconf=ind(2)
498           write (iout,*)
499           write (iout,*)
500           write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
501
502           etot=energy(0)
503           call enerprint(energy(0))
504           call briefout(it,etot)
505 c          if (minim) call briefout(it,etot)
506           if (refstr) then 
507             call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
508             write (istat,'(i5,19(f12.3))') iconf,
509      &     (energy(print_order(i)),i=1,nprint_ene),etot,
510      &     rms,frac,frac_nn,co
511           else
512             write (istat,'(i5,15(f12.3))') iconf,
513      &     (energy(print_order(i)),i=1,nprint_ene),etot
514           endif
515
516           imm=imm-1
517           if (read_cart) then
518             read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
519             call read_x(intin,*11)
520 #ifdef MPI
521 c Broadcast the order to compute internal coordinates to the slaves.
522             if (nfgtasks.gt.1)
523      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
524 #endif
525             call int_from_cart1(.false.)
526           else
527             read (intin,'(i5)',end=1101,err=1101) iconf
528             call read_angles(intin,*11)
529             call geom_to_var(nvar,varia)
530             call chainbuild
531           endif
532           n=n+1
533           imm=imm+1
534           ind(1)=1
535           ind(2)=n
536           ind(3)=0
537           ind(4)=0
538           ind(5)=0
539           ind(6)=0
540           call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
541      *                  ierr)
542           call mpi_send(varia,nvar,mpi_double_precision,man, 
543      *                  idreal,CG_COMM,ierr)
544           call mpi_send(ene0,1,mpi_double_precision,man,
545      *                  idreal,CG_COMM,ierr)
546           nf_mcmf=nf_mcmf+ind(4)
547           nmin=nmin+1
548         endif
549       enddo
550 11    continue
551       do j=1,imm
552         call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
553      *               CG_COMM,muster,ierr)
554         man=muster(mpi_source)
555         call mpi_recv(varia,nvar,mpi_double_precision, 
556      *               man,idreal,CG_COMM,muster,ierr)
557         call mpi_recv(ene,1,
558      *               mpi_double_precision,man,idreal,
559      *               CG_COMM,muster,ierr)
560         call mpi_recv(ene0,1,
561      *               mpi_double_precision,man,idreal,
562      *               CG_COMM,muster,ierr)
563
564         call var_to_geom(nvar,varia)
565         call chainbuild
566         call etotal(energy(0))
567         iconf=ind(2)
568         write (iout,*)
569         write (iout,*)
570         write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
571
572         etot=energy(0)
573         call enerprint(energy(0))
574         call briefout(it,etot)
575         if (refstr) then 
576           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
577           write (istat,'(i5,19(f12.3))') iconf,
578      &   (energy(print_order(i)),i=1,nprint_ene),etot,
579      &   rms,frac,frac_nn,co
580         else
581           write (istat,'(i5,15(f12.3))') iconf,
582      &    (energy(print_order(i)),i=1,nprint_ene),etot
583         endif
584         nmin=nmin+1
585       enddo
586 1101  continue
587       do i=1, nodes-1
588          ind(1)=0
589          ind(2)=0
590          ind(3)=0
591          ind(4)=0
592          ind(5)=0
593          ind(6)=0
594          call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
595      *                  ierr)
596       enddo
597 #else
598       close (intin)
599       open(intin,file=intinname,status='old')
600       write (istat,'(a5,20a12)')"#    ",
601      &   (wname(print_order(i)),i=1,nprint_ene)
602       write (istat,'("#    ",20(1pe12.4))')
603      &   (weights(print_order(i)),i=1,nprint_ene)
604       if (refstr) then
605         write (istat,'(a5,20a12)')"#    ",
606      &   (ename(print_order(i)),i=1,nprint_ene),
607      &   "ETOT total","RMSD","nat.contact","nnt.contact"
608       else
609         write (istat,'(a5,14a12)')"#    ",
610      &   (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
611       endif
612       do while (.not. eof)
613           if (read_cart) then
614             read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
615             call read_x(intin,*11)
616 #ifdef MPI
617 c Broadcast the order to compute internal coordinates to the slaves.
618             if (nfgtasks.gt.1)
619      &        call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
620 #endif
621             call int_from_cart1(.false.)
622           else
623             read (intin,'(i5)',end=1100,err=1100) iconf
624             call read_angles(intin,*11)
625             call geom_to_var(nvar,varia)
626             call chainbuild
627           endif
628         write (iout,'(a,i7)') 'Conformation #',iconf
629         if (minim) call minimize(etot,varia,iretcode,nfun)
630         call etotal(energy(0))
631
632         etot=energy(0)
633         call enerprint(energy(0))
634         if (minim) call briefout(it,etot) 
635         if (refstr) then 
636           call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
637           write (istat,'(i5,18(f12.3))') iconf,
638      &   (energy(print_order(i)),i=1,nprint_ene),
639      &   etot,rms,frac,frac_nn,co
640 cjlee end
641         else
642           write (istat,'(i5,14(f12.3))') iconf,
643      &   (energy(print_order(i)),i=1,nprint_ene),etot
644         endif
645       enddo
646    11 continue
647 #endif
648       return
649       end
650 c---------------------------------------------------------------------------
651       subroutine exec_checkgrad
652       implicit real*8 (a-h,o-z)
653       include 'DIMENSIONS'
654 #ifdef MPI
655       include 'mpif.h'
656 #endif
657       include 'COMMON.SETUP'
658       include 'COMMON.TIME1'
659       include 'COMMON.INTERACT'
660       include 'COMMON.NAMES'
661       include 'COMMON.GEO'
662       include 'COMMON.HEADER'
663       include 'COMMON.CONTROL'
664       include 'COMMON.CONTACTS'
665       include 'COMMON.CHAIN'
666       include 'COMMON.VAR'
667       include 'COMMON.IOUNITS'
668       include 'COMMON.FFIELD'
669       include 'COMMON.REMD'
670       include 'COMMON.MD'
671       include 'COMMON.SBRIDGE'
672       common /srutu/ icall
673       double precision energy(0:max_ene)
674 c      do i=2,nres
675 c        vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
676 c        if (itype(i).ne.10) 
677 c     &      vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
678 c      enddo
679       if (indpdb.eq.0) call chainbuild
680 c      do i=0,nres
681 c        do j=1,3
682 c          dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
683 c        enddo
684 c      enddo
685 c      do i=1,nres-1
686 c        if (itype(i).ne.10) then
687 c          do j=1,3
688 c            dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
689 c          enddo
690 c        endif
691 c      enddo
692 c      do j=1,3
693 c        dc(j,0)=ran_number(-0.2d0,0.2d0)
694 c      enddo
695       usampl=.true.
696       totT=1.d0
697       eq_time=0.0d0
698       call read_fragments
699       call chainbuild_cart
700       call cartprint
701       call intout
702       icall=1
703       call etotal(energy(0))
704       etot = energy(0)
705       call enerprint(energy(0))
706       write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
707       print *,'icheckgrad=',icheckgrad
708       goto (10,20,30) icheckgrad
709   10  call check_ecartint
710       return
711   20  call check_cartgrad
712       return
713   30  call check_eint
714       return
715       end
716 c---------------------------------------------------------------------------
717       subroutine exec_map
718 C Energy maps
719       call map_read
720       call map
721       return
722       end
723 c---------------------------------------------------------------------------
724       subroutine exec_CSA
725 #ifdef MPI
726       include "mpif.h"
727 #endif
728       include 'DIMENSIONS'
729       include 'COMMON.IOUNITS'
730 C Conformational Space Annealling programmed by Jooyoung Lee.
731 C This method works only with parallel machines!
732 #ifdef MPI
733       call together
734 #else
735       write (iout,*) "CSA works on parallel machines only"
736 #endif
737       return
738       end
739 c---------------------------------------------------------------------------
740       subroutine exec_softreg
741       include 'DIMENSIONS'
742       include 'COMMON.IOUNITS'
743       include 'COMMON.CONTROL'
744       double precision energy(0:max_ene)
745       call chainbuild
746       call etotal(energy(0))
747       call enerprint(energy(0))
748       if (.not.lsecondary) then
749         write(iout,*) 'Calling secondary structure recognition'
750         call secondary2(debug)
751       else
752         write(iout,*) 'Using secondary structure supplied in pdb'
753       endif
754
755       call softreg
756
757       call etotal(energy(0))
758       etot=energy(0)
759       call enerprint(energy(0))
760       call intout
761       call briefout(0,etot)
762       call secondary2(.true.)
763       if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
764       return
765       end