1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5 C Program to carry out conformational search of proteins in an united-residue C
8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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 c 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)'/
53 c call memmon_print_usage()
57 & write(iout,*)'### LAST MODIFIED 4/25/08 7:29PM by adam'
58 if (me.eq.king) call cinfo
59 C Read force field parameters and job setup data
62 c write (iout,*) "After readrtns"
65 if (me.eq.king .or. .not. out1file) then
67 & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
69 if (minim) write (iout,'(a)')
70 & 'Conformations will be energy-minimized.'
71 write (iout,'(80(1h*)/)')
75 if (modecalc.eq.-2) then
78 else if (modecalc.eq.-1) then
79 write(iout,*) "call check_sc_map next"
84 if (fg_rank.gt.0) then
85 C Fine-grain slaves just do energy and gradient components.
86 call ergastulum ! slave workhouse in Latin
89 if (modecalc.eq.0) then
90 write (iout,*) "Calling exec_eeval_or_minim"
92 call exec_eeval_or_minim
93 else if (modecalc.eq.1) then
95 else if (modecalc.eq.2) then
97 else if (modecalc.eq.3 .or. modecalc .eq.6) then
99 else if (modecalc.eq.4) then
100 call exec_mult_eeval_or_minim
101 else if (modecalc.eq.5) then
103 else if (ModeCalc.eq.7) then
105 else if (ModeCalc.eq.8) then
107 else if (modecalc.eq.11) then
109 else if (modecalc.eq.12) then
111 else if (modecalc.eq.14) then
115 write (iout,*) "Need a parallel version to run MREMD."
119 write (iout,'(a)') 'This calculation type is not supported',
125 if (fg_rank.eq.0) call finish_task
126 c call memmon_print_usage()
128 call print_detailed_timing
130 call MPI_Finalize(ierr)
133 call dajczas(tcpu(),hrtime,mintime,sectime)
134 stop '********** Program terminated normally.'
137 c--------------------------------------------------------------------------
144 include 'COMMON.SETUP'
145 include 'COMMON.CONTROL'
146 include 'COMMON.IOUNITS'
147 c if (me.eq.king .or. .not. out1file) then
148 c write (iout,*) "Calling chainbuild"
152 c if (me.eq.king .or. .not. out1file) then
153 c write (iout,*) "Calling MD"
159 c---------------------------------------------------------------------------
161 subroutine exec_MREMD
167 include 'COMMON.SETUP'
168 include 'COMMON.CONTROL'
169 include 'COMMON.IOUNITS'
170 include 'COMMON.REMD'
172 if (me.eq.king .or. .not. out1file)
173 & write (iout,*) "Calling chainbuild"
175 if (me.eq.king .or. .not. out1file)
176 & write (iout,*) "Calling REMD"
188 c---------------------------------------------------------------------------
189 subroutine exec_eeval_or_minim
195 include 'COMMON.SETUP'
196 include 'COMMON.TIME1'
197 include 'COMMON.INTERACT'
198 include 'COMMON.NAMES'
200 include 'COMMON.HEADER'
201 include 'COMMON.CONTROL'
202 c include 'COMMON.CONTACTS'
203 include 'COMMON.CHAIN'
205 include 'COMMON.IOUNITS'
206 include 'COMMON.FFIELD'
207 include 'COMMON.REMD'
209 include 'COMMON.SBRIDGE'
210 integer i,icall,iretcode,nfun
212 integer nharp,iharp(4,maxres/3)
215 double precision energy(0:n_ene),etot,etota
216 double precision energy_long(0:n_ene),energy_short(0:n_ene)
217 double precision rms,frac,frac_nn,co
218 double precision varia(maxvar)
219 double precision time00,time1,time_ene,evals
223 common /lbfgstat/ status,niter,nfun
226 if (indpdb.eq.0) call chainbuild
227 if (indpdb.ne.0) then
237 write (iout,*) "Energy evaluation/minimization"
239 c print *,'dc',dc(1,0),dc(2,0),dc(3,0)
241 print *,"Processor",myrank," after chainbuild"
243 call etotal_long(energy_long(0))
244 write (iout,*) "Printing long range energy"
245 call enerprint(energy_long(0))
246 call etotal_short(energy_short(0))
247 write (iout,*) "Printing short range energy"
248 call enerprint(energy_short(0))
250 energy(i)=energy_long(i)+energy_short(i)
251 c write (iout,*) i,energy_long(i),energy_short(i),energy(i)
253 write (iout,*) "Printing long+short range energy"
254 call enerprint(energy(0))
256 c write(iout,*)"before etotal"
258 call etotal(energy(0))
259 c write(iout,*)"after etotal"
262 time_ene=MPI_Wtime()-time00
264 time_ene=tcpu()-time00
266 write (iout,*) "Time for energy evaluation",time_ene
267 c print *,"after etotal"
270 call enerprint(energy(0))
271 call hairpin(.true.,nharp,iharp)
272 c print *,'after hairpin'
273 call secondary2(.true.)
274 c print *,'after secondary'
277 if (indpdb.ne.0 .and. .not.dccart) then
279 call chainbuild_extconf
280 call etotal(energy(0))
281 write (iout,*) "After bond regularization"
282 call enerprint(energy(0))
286 write (iout,*) 'Calling OVERLAP_SC'
287 call overlap_sc(fail)
288 write (iout,*) "After overlap_sc"
292 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
293 print *,'SC_move',nft_sc,etot
294 write(iout,*) 'SC_move',nft_sc,etot
298 print *, 'Calling MINIM_DC'
304 call minim_dc(etot,iretcode,nfun)
306 call geom_to_var(nvar,varia)
307 c print *,'Calling MINIMIZE.'
313 call minimize(etot,varia,iretcode,nfun)
316 print *,'LBFGS return code is',status,' eval ',nfun
318 print *,'SUMSL return code is',iretcode,' eval ',nfun
321 evals=nfun/(MPI_WTIME()-time1)
323 evals=nfun/(tcpu()-time1)
325 print *,'# eval/s',evals
326 print *,'refstr=',refstr
327 call hairpin(.false.,nharp,iharp)
328 print *,'after hairpin'
329 call secondary2(.true.)
330 print *,'after secondary'
331 call etotal(energy(0))
333 call enerprint(energy(0))
336 if (out_int) call briefout(0,etot)
338 cartname=prefix(:ilen(prefix))//'.x'
342 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
344 write (iout,'(a,a9)') 'LBFGS return code:',status
345 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
346 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
348 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
349 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
350 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
353 print *,'refstr=',refstr
354 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
355 call briefout(0,etot)
357 if (outpdb) call pdbout(etot,titel(:50),ipdb)
358 if (outmol2) call mol2out(etot,titel(:32))
361 c---------------------------------------------------------------------------
362 subroutine exec_regularize
368 include 'COMMON.SETUP'
369 include 'COMMON.TIME1'
370 include 'COMMON.INTERACT'
371 include 'COMMON.NAMES'
373 include 'COMMON.HEADER'
374 include 'COMMON.CONTROL'
375 c include 'COMMON.CONTACTS'
376 include 'COMMON.CHAIN'
378 include 'COMMON.IOUNITS'
379 include 'COMMON.FFIELD'
380 include 'COMMON.REMD'
382 include 'COMMON.SBRIDGE'
383 double precision energy(0:n_ene)
384 double precision etot,rms,frac,frac_nn,co
389 common /lbfgstat/ status,niter,nfun
395 call regularize(nct-nnt+1,etot,rms,cref(1,nnt),iretcode)
396 call etotal(energy(0))
397 energy(0)=energy(0)-energy(14)
399 call enerprint(energy(0))
401 call briefout(0,etot)
402 if (outpdb) call pdbout(etot,titel(:50),ipdb)
403 if (outmol2) call mol2out(etot,titel(:32))
404 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
406 write (iout,'(a,a9)') 'LBFGS return code:',status
408 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
412 c---------------------------------------------------------------------------
413 subroutine exec_thread
419 include "COMMON.SETUP"
423 c---------------------------------------------------------------------------
427 character*10 nodeinfo
429 double precision varia(maxvar)
433 include "COMMON.SETUP"
434 include 'COMMON.CONTROL'
438 if (modecalc.eq.3) then
444 if (modecalc.eq.3) then
455 c---------------------------------------------------------------------------
456 subroutine exec_mult_eeval_or_minim
461 integer muster(mpi_status_size)
464 include 'COMMON.SETUP'
465 include 'COMMON.TIME1'
466 include 'COMMON.INTERACT'
467 include 'COMMON.NAMES'
469 include 'COMMON.HEADER'
470 include 'COMMON.CONTROL'
471 c include 'COMMON.CONTACTS'
472 include 'COMMON.CHAIN'
474 include 'COMMON.IOUNITS'
475 include 'COMMON.FFIELD'
476 include 'COMMON.REMD'
478 include 'COMMON.SBRIDGE'
479 double precision varia(maxvar)
480 integer i,j,iconf,ind(6)
481 integer n,it,man,nf_mcmf,nmin,imm,mm,nft
482 double precision energy(0:max_ene),ene,etot,ene0
483 double precision rms,frac,frac_nn,co
484 double precision time
494 open(intin,file=intinname,status='old')
495 write (istat,'(a5,20a12)')"# ",
496 & (wname(print_order(i)),i=1,nprint_ene)
498 write (istat,'(a5,20a12)')"# ",
499 & (ename(print_order(i)),i=1,nprint_ene),
500 & "ETOT total","RMSD","nat.contact","nnt.contact"
502 write (istat,'(a5,20a12)')"# ",
503 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
509 read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
510 call read_x(intin,*11)
512 c Broadcast the order to compute internal coordinates to the slaves.
514 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
516 call int_from_cart1(.false.)
518 read (intin,'(i5)',end=1100,err=1100) iconf
519 call read_angles(intin,*11)
520 call geom_to_var(nvar,varia)
523 write (iout,'(a,i7)') 'Conformation #',iconf
524 call etotal(energy(0))
525 call briefout(iconf,energy(0))
526 call enerprint(energy(0))
529 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
530 write (istat,'(i5,20(f12.3))') iconf,
531 & (energy(print_order(i)),i=1,nprint_ene),etot,
532 & rms,frac,frac_nn,co
535 write (istat,'(i5,16(f12.3))') iconf,
536 & (energy(print_order(i)),i=1,nprint_ene),etot
552 if (mm.lt.nodes) then
554 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
555 call read_x(intin,*11)
557 c Broadcast the order to compute internal coordinates to the slaves.
559 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
561 call int_from_cart1(.false.)
563 read (intin,'(i5)',end=11,err=11) iconf
564 call read_angles(intin,*11)
565 call geom_to_var(nvar,varia)
568 write (iout,'(a,i7)') 'Conformation #',iconf
578 call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
580 call mpi_send(varia,nvar,mpi_double_precision,mm,
581 * idreal,CG_COMM,ierr)
582 call mpi_send(ene0,1,mpi_double_precision,mm,
583 * idreal,CG_COMM,ierr)
584 c print *,'task ',n,' sent to worker ',mm,nvar
586 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
587 * CG_COMM,muster,ierr)
588 man=muster(mpi_source)
589 c print *,'receiving result from worker ',man,' (',iii1,iii,')'
590 call mpi_recv(varia,nvar,mpi_double_precision,
591 * man,idreal,CG_COMM,muster,ierr)
593 * mpi_double_precision,man,idreal,
594 * CG_COMM,muster,ierr)
595 call mpi_recv(ene0,1,
596 * mpi_double_precision,man,idreal,
597 * CG_COMM,muster,ierr)
598 c print *,'result received from worker ',man,' sending now'
600 call var_to_geom(nvar,varia)
602 call etotal(energy(0))
606 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
609 call enerprint(energy(0))
610 call briefout(it,etot)
611 c if (minim) call briefout(it,etot)
613 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
614 write (istat,'(i5,19(f12.3))') iconf,
615 & (energy(print_order(i)),i=1,nprint_ene),etot,
616 & rms,frac,frac_nn,co
618 write (istat,'(i5,15(f12.3))') iconf,
619 & (energy(print_order(i)),i=1,nprint_ene),etot
624 read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
625 call read_x(intin,*11)
627 c Broadcast the order to compute internal coordinates to the slaves.
629 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
631 call int_from_cart1(.false.)
633 read (intin,'(i5)',end=1101,err=1101) iconf
634 call read_angles(intin,*11)
635 call geom_to_var(nvar,varia)
646 call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
648 call mpi_send(varia,nvar,mpi_double_precision,man,
649 * idreal,CG_COMM,ierr)
650 call mpi_send(ene0,1,mpi_double_precision,man,
651 * idreal,CG_COMM,ierr)
652 nf_mcmf=nf_mcmf+ind(4)
658 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
659 * CG_COMM,muster,ierr)
660 man=muster(mpi_source)
661 call mpi_recv(varia,nvar,mpi_double_precision,
662 * man,idreal,CG_COMM,muster,ierr)
664 * mpi_double_precision,man,idreal,
665 * CG_COMM,muster,ierr)
666 call mpi_recv(ene0,1,
667 * mpi_double_precision,man,idreal,
668 * CG_COMM,muster,ierr)
670 call var_to_geom(nvar,varia)
672 call etotal(energy(0))
676 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
679 call enerprint(energy(0))
680 call briefout(it,etot)
682 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
683 write (istat,'(i5,19(f12.3))') iconf,
684 & (energy(print_order(i)),i=1,nprint_ene),etot,
685 & rms,frac,frac_nn,co
687 write (istat,'(i5,15(f12.3))') iconf,
688 & (energy(print_order(i)),i=1,nprint_ene),etot
700 call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
705 open(intin,file=intinname,status='old')
706 write (istat,'(a5,20a12)')"# ",
707 & (wname(print_order(i)),i=1,nprint_ene)
708 write (istat,'("# ",20(1pe12.4))')
709 & (weights(print_order(i)),i=1,nprint_ene)
711 write (istat,'(a5,20a12)')"# ",
712 & (ename(print_order(i)),i=1,nprint_ene),
713 & "ETOT total","RMSD","nat.contact","nnt.contact"
715 write (istat,'(a5,14a12)')"# ",
716 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
720 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
721 call read_x(intin,*11)
723 c Broadcast the order to compute internal coordinates to the slaves.
725 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
727 call int_from_cart1(.false.)
729 read (intin,'(i5)',end=11,err=11) iconf
730 call read_angles(intin,*11)
731 call geom_to_var(nvar,varia)
734 write (iout,'(a,i7)') 'Conformation #',iconf
735 if (minim) call minimize(etot,varia,iretcode,nfun)
736 call etotal(energy(0))
739 call enerprint(energy(0))
740 if (minim) call briefout(it,etot)
742 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
743 write (istat,'(i5,18(f12.3))') iconf,
744 & (energy(print_order(i)),i=1,nprint_ene),
745 & etot,rms,frac,frac_nn,co
748 write (istat,'(i5,14(f12.3))') iconf,
749 & (energy(print_order(i)),i=1,nprint_ene),etot
756 c---------------------------------------------------------------------------
757 subroutine exec_checkgrad
763 include 'COMMON.SETUP'
764 include 'COMMON.TIME1'
765 include 'COMMON.INTERACT'
766 include 'COMMON.NAMES'
768 include 'COMMON.HEADER'
769 include 'COMMON.CONTROL'
770 c include 'COMMON.CONTACTS'
771 include 'COMMON.CHAIN'
773 include 'COMMON.IOUNITS'
774 include 'COMMON.FFIELD'
775 include 'COMMON.REMD'
777 include 'COMMON.QRESTR'
778 include 'COMMON.SBRIDGE'
781 double precision energy(0:max_ene)
784 c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
785 c if (itype(i).ne.10)
786 c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
788 if (indpdb.eq.0) call chainbuild
791 c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
795 c if (itype(i).ne.10) then
797 c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
802 c dc(j,0)=ran_number(-0.2d0,0.2d0)
815 print *, "AFTER read fragments"
816 write (iout,*) "iset",iset
818 write(iout,*) "fragment, weights, q0:"
820 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
821 & ifrag_back(2,i,iset),
822 & wfrag_back(1,i,iset),qin_back(1,i,iset),
823 & wfrag_back(2,i,iset),qin_back(2,i,iset),
824 & wfrag_back(3,i,iset),qin_back(3,i,iset)
827 write(iout,*) "fragment, weights:"
829 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
830 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
831 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
838 call rescale_weights(t_bath)
840 print *,"chainbuild_cart"
842 print *,"After cartprint"
845 print *,"before ETOT"
846 write (iout,*) "usampl",usampl
847 call etotal(energy(0))
849 call enerprint(energy(0))
850 write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
851 print *,'icheckgrad=',icheckgrad
853 goto (10,20,30) icheckgrad
854 10 call check_ecartint
857 & "Checking the gradient of Cartesian coordinates disabled."
862 c---------------------------------------------------------------------------
869 c---------------------------------------------------------------------------
876 include 'COMMON.IOUNITS'
877 C Conformational Space Annealling programmed by Jooyoung Lee.
878 C This method works only with parallel machines!
882 write (iout,*) "CSA works on parallel machines only"
886 c---------------------------------------------------------------------------
887 subroutine exec_softreg
890 include 'COMMON.IOUNITS'
891 include 'COMMON.CONTROL'
892 double precision energy(0:max_ene),etot
893 double precision rms,frac,frac_nn,co
895 call etotal(energy(0))
896 call enerprint(energy(0))
897 if (.not.lsecondary) then
898 write(iout,*) 'Calling secondary structure recognition'
899 call secondary2(.true.)
901 write(iout,*) 'Using secondary structure supplied in pdb'
906 call etotal(energy(0))
908 call enerprint(energy(0))
910 call briefout(0,etot)
911 call secondary2(.true.)
912 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)