1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5 C Program to carry out conformational search of proteins in an united-residue C
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
9 implicit real*8 (a-h,o-z)
15 include 'COMMON.SETUP'
17 include 'COMMON.TIME1'
18 include 'COMMON.INTERACT'
19 include 'COMMON.NAMES'
21 include 'COMMON.HEADER'
22 include 'COMMON.CONTROL'
23 include 'COMMON.CONTACTS'
24 include 'COMMON.CHAIN'
26 include 'COMMON.IOUNITS'
27 include 'COMMON.FFIELD'
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)',
45 & 'Soft regularization of PDB structure',
46 & 'Mesoscopic molecular dynamics (MD) ',
48 & 'Replica exchange molecular dynamics (REMD)'/
51 c call memmon_print_usage()
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
60 c write (iout,*) "After readrtns"
63 if (me.eq.king .or. .not. out1file) then
65 & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
67 if (minim) write (iout,'(a)')
68 & 'Conformations will be energy-minimized.'
69 write (iout,'(80(1h*)/)')
73 if (modecalc.eq.-2) then
76 else if (modecalc.eq.-1) then
77 write(iout,*) "call check_sc_map next"
82 if (fg_rank.gt.0) then
83 C Fine-grain slaves just do energy and gradient components.
84 call ergastulum ! slave workhouse in Latin
87 if (modecalc.eq.0) then
88 write (iout,*) "Calling exec_eeval_or_minim"
90 call exec_eeval_or_minim
91 else if (modecalc.eq.1) then
93 else if (modecalc.eq.2) then
95 else if (modecalc.eq.3 .or. modecalc .eq.6) then
97 else if (modecalc.eq.4) then
98 call exec_mult_eeval_or_minim
99 else if (modecalc.eq.5) then
101 else if (ModeCalc.eq.7) then
103 else if (ModeCalc.eq.8) then
105 else if (modecalc.eq.11) then
107 else if (modecalc.eq.12) then
109 else if (modecalc.eq.14) then
113 write (iout,*) "Need a parallel version to run MREMD."
117 write (iout,'(a)') 'This calculation type is not supported',
123 if (fg_rank.eq.0) call finish_task
124 c call memmon_print_usage()
126 call print_detailed_timing
128 call MPI_Finalize(ierr)
131 call dajczas(tcpu(),hrtime,mintime,sectime)
132 stop '********** Program terminated normally.'
135 c--------------------------------------------------------------------------
141 include 'COMMON.SETUP'
142 include 'COMMON.CONTROL'
143 include 'COMMON.IOUNITS'
144 c if (me.eq.king .or. .not. out1file) then
145 c write (iout,*) "Calling chainbuild"
149 c if (me.eq.king .or. .not. out1file) then
150 c write (iout,*) "Calling MD"
156 c---------------------------------------------------------------------------
158 subroutine exec_MREMD
163 include 'COMMON.SETUP'
164 include 'COMMON.CONTROL'
165 include 'COMMON.IOUNITS'
166 include 'COMMON.REMD'
167 if (me.eq.king .or. .not. out1file)
168 & write (iout,*) "Calling chainbuild"
170 if (me.eq.king .or. .not. out1file)
171 & write (iout,*) "Calling REMD"
183 c---------------------------------------------------------------------------
184 subroutine exec_eeval_or_minim
185 implicit real*8 (a-h,o-z)
190 include 'COMMON.SETUP'
191 include 'COMMON.TIME1'
192 include 'COMMON.INTERACT'
193 include 'COMMON.NAMES'
195 include 'COMMON.HEADER'
196 include 'COMMON.CONTROL'
197 include 'COMMON.CONTACTS'
198 include 'COMMON.CHAIN'
200 include 'COMMON.IOUNITS'
201 include 'COMMON.FFIELD'
202 include 'COMMON.REMD'
204 include 'COMMON.SBRIDGE'
206 double precision energy(0:n_ene)
207 double precision energy_long(0:n_ene),energy_short(0:n_ene)
208 double precision varia(maxvar)
209 if (indpdb.eq.0) call chainbuild
210 if (indpdb.ne.0) then
220 write (iout,*) "Energy evaluation/minimization"
222 c print *,'dc',dc(1,0),dc(2,0),dc(3,0)
224 print *,"Processor",myrank," after chainbuild"
226 call etotal_long(energy_long(0))
227 write (iout,*) "Printing long range energy"
228 call enerprint(energy_long(0))
229 call etotal_short(energy_short(0))
230 write (iout,*) "Printing short range energy"
231 call enerprint(energy_short(0))
233 energy(i)=energy_long(i)+energy_short(i)
234 c write (iout,*) i,energy_long(i),energy_short(i),energy(i)
236 write (iout,*) "Printing long+short range energy"
237 call enerprint(energy(0))
239 c write(iout,*)"before etotal"
241 call etotal(energy(0))
242 c write(iout,*)"after etotal"
245 time_ene=MPI_Wtime()-time00
247 time_ene=tcpu()-time00
249 write (iout,*) "Time for energy evaluation",time_ene
250 print *,"after etotal"
253 call enerprint(energy(0))
254 call hairpin(.true.,nharp,iharp)
255 print *,'after hairpin'
256 call secondary2(.true.)
257 print *,'after secondary'
261 print *, 'Calling OVERLAP_SC'
262 call overlap_sc(fail)
263 print *,"After overlap_sc"
267 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
268 print *,'SC_move',nft_sc,etot
269 write(iout,*) 'SC_move',nft_sc,etot
273 print *, 'Calling MINIM_DC'
279 call minim_dc(etot,iretcode,nfun)
281 if (indpdb.ne.0) then
283 call chainbuild_extconf
285 call geom_to_var(nvar,varia)
286 print *,'Calling MINIMIZE.'
292 call minimize(etot,varia,iretcode,nfun)
294 print *,'SUMSL return code is',iretcode,' eval ',nfun
296 evals=nfun/(MPI_WTIME()-time1)
298 evals=nfun/(tcpu()-time1)
300 print *,'# eval/s',evals
301 print *,'refstr=',refstr
302 call hairpin(.false.,nharp,iharp)
303 print *,'after hairpin'
304 call secondary2(.true.)
305 print *,'after secondary'
306 call etotal(energy(0))
308 call enerprint(energy(0))
311 call briefout(0,etot)
312 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
313 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
314 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
315 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
317 print *,'refstr=',refstr
318 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
319 call briefout(0,etot)
321 if (outpdb) call pdbout(etot,titel(:50),ipdb)
322 if (outmol2) call mol2out(etot,titel(:32))
325 c---------------------------------------------------------------------------
326 subroutine exec_regularize
327 implicit real*8 (a-h,o-z)
332 include 'COMMON.SETUP'
333 include 'COMMON.TIME1'
334 include 'COMMON.INTERACT'
335 include 'COMMON.NAMES'
337 include 'COMMON.HEADER'
338 include 'COMMON.CONTROL'
339 include 'COMMON.CONTACTS'
340 include 'COMMON.CHAIN'
342 include 'COMMON.IOUNITS'
343 include 'COMMON.FFIELD'
344 include 'COMMON.REMD'
346 include 'COMMON.SBRIDGE'
347 double precision energy(0:n_ene)
352 call regularize(nct-nnt+1,etot,rms,cref(1,nnt),iretcode)
353 call etotal(energy(0))
354 energy(0)=energy(0)-energy(14)
356 call enerprint(energy(0))
358 call briefout(0,etot)
359 if (outpdb) call pdbout(etot,titel(:50),ipdb)
360 if (outmol2) call mol2out(etot,titel(:32))
361 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
362 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
365 c---------------------------------------------------------------------------
366 subroutine exec_thread
371 include "COMMON.SETUP"
375 c---------------------------------------------------------------------------
377 implicit real*8 (a-h,o-z)
379 character*10 nodeinfo
380 double precision varia(maxvar)
384 include "COMMON.SETUP"
385 include 'COMMON.CONTROL'
389 if (modecalc.eq.3) then
395 if (modecalc.eq.3) then
406 c---------------------------------------------------------------------------
407 subroutine exec_mult_eeval_or_minim
408 implicit real*8 (a-h,o-z)
412 dimension muster(mpi_status_size)
414 include 'COMMON.SETUP'
415 include 'COMMON.TIME1'
416 include 'COMMON.INTERACT'
417 include 'COMMON.NAMES'
419 include 'COMMON.HEADER'
420 include 'COMMON.CONTROL'
421 include 'COMMON.CONTACTS'
422 include 'COMMON.CHAIN'
424 include 'COMMON.IOUNITS'
425 include 'COMMON.FFIELD'
426 include 'COMMON.REMD'
428 include 'COMMON.SBRIDGE'
429 double precision varia(maxvar)
431 double precision energy(0:max_ene)
441 open(intin,file=intinname,status='old')
442 write (istat,'(a5,20a12)')"# ",
443 & (wname(print_order(i)),i=1,nprint_ene)
445 write (istat,'(a5,20a12)')"# ",
446 & (ename(print_order(i)),i=1,nprint_ene),
447 & "ETOT total","RMSD","nat.contact","nnt.contact"
449 write (istat,'(a5,20a12)')"# ",
450 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
456 read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
457 call read_x(intin,*11)
459 c Broadcast the order to compute internal coordinates to the slaves.
461 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
463 call int_from_cart1(.false.)
465 read (intin,'(i5)',end=1100,err=1100) iconf
466 call read_angles(intin,*11)
467 call geom_to_var(nvar,varia)
470 write (iout,'(a,i7)') 'Conformation #',iconf
471 call etotal(energy(0))
472 call briefout(iconf,energy(0))
473 call enerprint(energy(0))
476 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
477 write (istat,'(i5,20(f12.3))') iconf,
478 & (energy(print_order(i)),i=1,nprint_ene),etot,
479 & rms,frac,frac_nn,co
482 write (istat,'(i5,16(f12.3))') iconf,
483 & (energy(print_order(i)),i=1,nprint_ene),etot
499 if (mm.lt.nodes) then
501 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
502 call read_x(intin,*11)
504 c Broadcast the order to compute internal coordinates to the slaves.
506 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
508 call int_from_cart1(.false.)
510 read (intin,'(i5)',end=11,err=11) iconf
511 call read_angles(intin,*11)
512 call geom_to_var(nvar,varia)
515 write (iout,'(a,i7)') 'Conformation #',iconf
525 call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
527 call mpi_send(varia,nvar,mpi_double_precision,mm,
528 * idreal,CG_COMM,ierr)
529 call mpi_send(ene0,1,mpi_double_precision,mm,
530 * idreal,CG_COMM,ierr)
531 c print *,'task ',n,' sent to worker ',mm,nvar
533 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
534 * CG_COMM,muster,ierr)
535 man=muster(mpi_source)
536 c print *,'receiving result from worker ',man,' (',iii1,iii,')'
537 call mpi_recv(varia,nvar,mpi_double_precision,
538 * man,idreal,CG_COMM,muster,ierr)
540 * mpi_double_precision,man,idreal,
541 * CG_COMM,muster,ierr)
542 call mpi_recv(ene0,1,
543 * mpi_double_precision,man,idreal,
544 * CG_COMM,muster,ierr)
545 c print *,'result received from worker ',man,' sending now'
547 call var_to_geom(nvar,varia)
549 call etotal(energy(0))
553 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
556 call enerprint(energy(0))
557 call briefout(it,etot)
558 c if (minim) call briefout(it,etot)
560 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
561 write (istat,'(i5,19(f12.3))') iconf,
562 & (energy(print_order(i)),i=1,nprint_ene),etot,
563 & rms,frac,frac_nn,co
565 write (istat,'(i5,15(f12.3))') iconf,
566 & (energy(print_order(i)),i=1,nprint_ene),etot
571 read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
572 call read_x(intin,*11)
574 c Broadcast the order to compute internal coordinates to the slaves.
576 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
578 call int_from_cart1(.false.)
580 read (intin,'(i5)',end=1101,err=1101) iconf
581 call read_angles(intin,*11)
582 call geom_to_var(nvar,varia)
593 call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
595 call mpi_send(varia,nvar,mpi_double_precision,man,
596 * idreal,CG_COMM,ierr)
597 call mpi_send(ene0,1,mpi_double_precision,man,
598 * idreal,CG_COMM,ierr)
599 nf_mcmf=nf_mcmf+ind(4)
605 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
606 * CG_COMM,muster,ierr)
607 man=muster(mpi_source)
608 call mpi_recv(varia,nvar,mpi_double_precision,
609 * man,idreal,CG_COMM,muster,ierr)
611 * mpi_double_precision,man,idreal,
612 * CG_COMM,muster,ierr)
613 call mpi_recv(ene0,1,
614 * mpi_double_precision,man,idreal,
615 * CG_COMM,muster,ierr)
617 call var_to_geom(nvar,varia)
619 call etotal(energy(0))
623 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
626 call enerprint(energy(0))
627 call briefout(it,etot)
629 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
630 write (istat,'(i5,19(f12.3))') iconf,
631 & (energy(print_order(i)),i=1,nprint_ene),etot,
632 & rms,frac,frac_nn,co
634 write (istat,'(i5,15(f12.3))') iconf,
635 & (energy(print_order(i)),i=1,nprint_ene),etot
647 call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
652 open(intin,file=intinname,status='old')
653 write (istat,'(a5,20a12)')"# ",
654 & (wname(print_order(i)),i=1,nprint_ene)
655 write (istat,'("# ",20(1pe12.4))')
656 & (weights(print_order(i)),i=1,nprint_ene)
658 write (istat,'(a5,20a12)')"# ",
659 & (ename(print_order(i)),i=1,nprint_ene),
660 & "ETOT total","RMSD","nat.contact","nnt.contact"
662 write (istat,'(a5,14a12)')"# ",
663 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
667 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
668 call read_x(intin,*11)
670 c Broadcast the order to compute internal coordinates to the slaves.
672 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
674 call int_from_cart1(.false.)
676 read (intin,'(i5)',end=11,err=11) iconf
677 call read_angles(intin,*11)
678 call geom_to_var(nvar,varia)
681 write (iout,'(a,i7)') 'Conformation #',iconf
682 if (minim) call minimize(etot,varia,iretcode,nfun)
683 call etotal(energy(0))
686 call enerprint(energy(0))
687 if (minim) call briefout(it,etot)
689 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
690 write (istat,'(i5,18(f12.3))') iconf,
691 & (energy(print_order(i)),i=1,nprint_ene),
692 & etot,rms,frac,frac_nn,co
695 write (istat,'(i5,14(f12.3))') iconf,
696 & (energy(print_order(i)),i=1,nprint_ene),etot
703 c---------------------------------------------------------------------------
704 subroutine exec_checkgrad
705 implicit real*8 (a-h,o-z)
710 include 'COMMON.SETUP'
711 include 'COMMON.TIME1'
712 include 'COMMON.INTERACT'
713 include 'COMMON.NAMES'
715 include 'COMMON.HEADER'
716 include 'COMMON.CONTROL'
717 include 'COMMON.CONTACTS'
718 include 'COMMON.CHAIN'
720 include 'COMMON.IOUNITS'
721 include 'COMMON.FFIELD'
722 include 'COMMON.REMD'
724 include 'COMMON.SBRIDGE'
726 double precision energy(0:max_ene)
729 c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
730 c if (itype(i).ne.10)
731 c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
733 if (indpdb.eq.0) call chainbuild
736 c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
740 c if (itype(i).ne.10) then
742 c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
747 c dc(j,0)=ran_number(-0.2d0,0.2d0)
760 print *, "AFTER read fragments"
761 write (iout,*) "iset",iset
763 write(iout,*) "fragment, weights, q0:"
765 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
766 & ifrag_back(2,i,iset),
767 & wfrag_back(1,i,iset),qin_back(1,i,iset),
768 & wfrag_back(2,i,iset),qin_back(2,i,iset),
769 & wfrag_back(3,i,iset),qin_back(3,i,iset)
772 write(iout,*) "fragment, weights:"
774 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
775 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
776 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
783 call rescale_weights(t_bath)
785 print *,"chainbuild_cart"
787 print *,"After cartprint"
790 print *,"before ETOT"
791 write (iout,*) "usampl",usampl
792 call etotal(energy(0))
794 call enerprint(energy(0))
795 write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
796 print *,'icheckgrad=',icheckgrad
798 goto (10,20,30) icheckgrad
799 10 call check_ecartint
801 20 call check_cartgrad
806 c---------------------------------------------------------------------------
813 c---------------------------------------------------------------------------
819 include 'COMMON.IOUNITS'
820 C Conformational Space Annealling programmed by Jooyoung Lee.
821 C This method works only with parallel machines!
825 write (iout,*) "CSA works on parallel machines only"
829 c---------------------------------------------------------------------------
830 subroutine exec_softreg
832 include 'COMMON.IOUNITS'
833 include 'COMMON.CONTROL'
834 double precision energy(0:max_ene)
836 call etotal(energy(0))
837 call enerprint(energy(0))
838 if (.not.lsecondary) then
839 write(iout,*) 'Calling secondary structure recognition'
840 call secondary2(debug)
842 write(iout,*) 'Using secondary structure supplied in pdb'
847 call etotal(energy(0))
849 call enerprint(energy(0))
851 call briefout(0,etot)
852 call secondary2(.true.)
853 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)