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 (outpdb) call pdbout(etot,titel(:50),ipdb)
343 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
345 write (iout,'(a,a9)') 'LBFGS return code:',status
346 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
347 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
349 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
350 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
351 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
354 print *,'refstr=',refstr
355 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
356 call briefout(0,etot)
358 if (outpdb) call pdbout(etot,titel(:50),ipdb)
359 if (outmol2) call mol2out(etot,titel(:32))
362 c---------------------------------------------------------------------------
363 subroutine exec_regularize
369 include 'COMMON.SETUP'
370 include 'COMMON.TIME1'
371 include 'COMMON.INTERACT'
372 include 'COMMON.NAMES'
374 include 'COMMON.HEADER'
375 include 'COMMON.CONTROL'
376 c include 'COMMON.CONTACTS'
377 include 'COMMON.CHAIN'
379 include 'COMMON.IOUNITS'
380 include 'COMMON.FFIELD'
381 include 'COMMON.REMD'
383 include 'COMMON.SBRIDGE'
384 double precision energy(0:n_ene)
385 double precision etot,rms,frac,frac_nn,co
390 common /lbfgstat/ status,niter,nfun
396 call regularize(nct-nnt+1,etot,rms,cref(1,nnt),iretcode)
397 call etotal(energy(0))
398 energy(0)=energy(0)-energy(14)
400 call enerprint(energy(0))
402 call briefout(0,etot)
403 if (outpdb) call pdbout(etot,titel(:50),ipdb)
404 if (outmol2) call mol2out(etot,titel(:32))
405 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
407 write (iout,'(a,a9)') 'LBFGS return code:',status
409 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
413 c---------------------------------------------------------------------------
414 subroutine exec_thread
420 include "COMMON.SETUP"
424 c---------------------------------------------------------------------------
428 character*10 nodeinfo
430 double precision varia(maxvar)
434 include "COMMON.SETUP"
435 include 'COMMON.CONTROL'
439 if (modecalc.eq.3) then
445 if (modecalc.eq.3) then
456 c---------------------------------------------------------------------------
457 subroutine exec_mult_eeval_or_minim
462 integer muster(mpi_status_size)
465 include 'COMMON.SETUP'
466 include 'COMMON.TIME1'
467 include 'COMMON.INTERACT'
468 include 'COMMON.NAMES'
470 include 'COMMON.HEADER'
471 include 'COMMON.CONTROL'
472 c include 'COMMON.CONTACTS'
473 include 'COMMON.CHAIN'
475 include 'COMMON.IOUNITS'
476 include 'COMMON.FFIELD'
477 include 'COMMON.REMD'
479 include 'COMMON.SBRIDGE'
480 double precision varia(maxvar)
481 integer i,j,iconf,ind(6)
482 integer n,it,man,nf_mcmf,nmin,imm,mm,nft
483 double precision energy(0:max_ene),ene,etot,ene0
484 double precision rms,frac,frac_nn,co
485 double precision time
495 open(intin,file=intinname,status='old')
496 write (istat,'(a5,20a12)')"# ",
497 & (wname(print_order(i)),i=1,nprint_ene)
499 write (istat,'(a5,20a12)')"# ",
500 & (ename(print_order(i)),i=1,nprint_ene),
501 & "ETOT total","RMSD","nat.contact","nnt.contact"
503 write (istat,'(a5,20a12)')"# ",
504 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
510 read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
511 call read_x(intin,*11)
513 c Broadcast the order to compute internal coordinates to the slaves.
515 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
517 call int_from_cart1(.false.)
519 read (intin,'(i5)',end=1100,err=1100) iconf
520 call read_angles(intin,*11)
521 call geom_to_var(nvar,varia)
524 write (iout,'(a,i7)') 'Conformation #',iconf
525 call etotal(energy(0))
526 call briefout(iconf,energy(0))
527 call enerprint(energy(0))
530 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
531 write (istat,'(i5,20(f12.3))') iconf,
532 & (energy(print_order(i)),i=1,nprint_ene),etot,
533 & rms,frac,frac_nn,co
536 write (istat,'(i5,16(f12.3))') iconf,
537 & (energy(print_order(i)),i=1,nprint_ene),etot
553 if (mm.lt.nodes) then
555 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
556 call read_x(intin,*11)
558 c Broadcast the order to compute internal coordinates to the slaves.
560 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
562 call int_from_cart1(.false.)
564 read (intin,'(i5)',end=11,err=11) iconf
565 call read_angles(intin,*11)
566 call geom_to_var(nvar,varia)
569 write (iout,'(a,i7)') 'Conformation #',iconf
579 call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
581 call mpi_send(varia,nvar,mpi_double_precision,mm,
582 * idreal,CG_COMM,ierr)
583 call mpi_send(ene0,1,mpi_double_precision,mm,
584 * idreal,CG_COMM,ierr)
585 c print *,'task ',n,' sent to worker ',mm,nvar
587 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
588 * CG_COMM,muster,ierr)
589 man=muster(mpi_source)
590 c print *,'receiving result from worker ',man,' (',iii1,iii,')'
591 call mpi_recv(varia,nvar,mpi_double_precision,
592 * man,idreal,CG_COMM,muster,ierr)
594 * mpi_double_precision,man,idreal,
595 * CG_COMM,muster,ierr)
596 call mpi_recv(ene0,1,
597 * mpi_double_precision,man,idreal,
598 * CG_COMM,muster,ierr)
599 c print *,'result received from worker ',man,' sending now'
601 call var_to_geom(nvar,varia)
603 call etotal(energy(0))
607 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
610 call enerprint(energy(0))
611 call briefout(it,etot)
612 c if (minim) call briefout(it,etot)
614 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
615 write (istat,'(i5,19(f12.3))') iconf,
616 & (energy(print_order(i)),i=1,nprint_ene),etot,
617 & rms,frac,frac_nn,co
619 write (istat,'(i5,15(f12.3))') iconf,
620 & (energy(print_order(i)),i=1,nprint_ene),etot
625 read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
626 call read_x(intin,*11)
628 c Broadcast the order to compute internal coordinates to the slaves.
630 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
632 call int_from_cart1(.false.)
634 read (intin,'(i5)',end=1101,err=1101) iconf
635 call read_angles(intin,*11)
636 call geom_to_var(nvar,varia)
647 call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
649 call mpi_send(varia,nvar,mpi_double_precision,man,
650 * idreal,CG_COMM,ierr)
651 call mpi_send(ene0,1,mpi_double_precision,man,
652 * idreal,CG_COMM,ierr)
653 nf_mcmf=nf_mcmf+ind(4)
659 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
660 * CG_COMM,muster,ierr)
661 man=muster(mpi_source)
662 call mpi_recv(varia,nvar,mpi_double_precision,
663 * man,idreal,CG_COMM,muster,ierr)
665 * mpi_double_precision,man,idreal,
666 * CG_COMM,muster,ierr)
667 call mpi_recv(ene0,1,
668 * mpi_double_precision,man,idreal,
669 * CG_COMM,muster,ierr)
671 call var_to_geom(nvar,varia)
673 call etotal(energy(0))
677 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
680 call enerprint(energy(0))
681 call briefout(it,etot)
683 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
684 write (istat,'(i5,19(f12.3))') iconf,
685 & (energy(print_order(i)),i=1,nprint_ene),etot,
686 & rms,frac,frac_nn,co
688 write (istat,'(i5,15(f12.3))') iconf,
689 & (energy(print_order(i)),i=1,nprint_ene),etot
701 call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
706 open(intin,file=intinname,status='old')
707 write (istat,'(a5,20a12)')"# ",
708 & (wname(print_order(i)),i=1,nprint_ene)
709 write (istat,'("# ",20(1pe12.4))')
710 & (weights(print_order(i)),i=1,nprint_ene)
712 write (istat,'(a5,20a12)')"# ",
713 & (ename(print_order(i)),i=1,nprint_ene),
714 & "ETOT total","RMSD","nat.contact","nnt.contact"
716 write (istat,'(a5,14a12)')"# ",
717 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
721 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
722 call read_x(intin,*11)
724 c Broadcast the order to compute internal coordinates to the slaves.
726 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
728 call int_from_cart1(.false.)
730 read (intin,'(i5)',end=11,err=11) iconf
731 call read_angles(intin,*11)
732 call geom_to_var(nvar,varia)
735 write (iout,'(a,i7)') 'Conformation #',iconf
736 if (minim) call minimize(etot,varia,iretcode,nfun)
737 call etotal(energy(0))
740 call enerprint(energy(0))
741 if (minim) call briefout(it,etot)
743 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
744 write (istat,'(i5,18(f12.3))') iconf,
745 & (energy(print_order(i)),i=1,nprint_ene),
746 & etot,rms,frac,frac_nn,co
749 write (istat,'(i5,14(f12.3))') iconf,
750 & (energy(print_order(i)),i=1,nprint_ene),etot
757 c---------------------------------------------------------------------------
758 subroutine exec_checkgrad
764 include 'COMMON.SETUP'
765 include 'COMMON.TIME1'
766 include 'COMMON.INTERACT'
767 include 'COMMON.NAMES'
769 include 'COMMON.HEADER'
770 include 'COMMON.CONTROL'
771 c include 'COMMON.CONTACTS'
772 include 'COMMON.CHAIN'
774 include 'COMMON.IOUNITS'
775 include 'COMMON.FFIELD'
776 include 'COMMON.REMD'
778 include 'COMMON.QRESTR'
779 include 'COMMON.SBRIDGE'
782 double precision energy(0:max_ene)
785 c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
786 c if (itype(i).ne.10)
787 c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
789 if (indpdb.eq.0) call chainbuild
792 c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
796 c if (itype(i).ne.10) then
798 c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
803 c dc(j,0)=ran_number(-0.2d0,0.2d0)
816 print *, "AFTER read fragments"
817 write (iout,*) "iset",iset
819 write(iout,*) "fragment, weights, q0:"
821 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
822 & ifrag_back(2,i,iset),
823 & wfrag_back(1,i,iset),qin_back(1,i,iset),
824 & wfrag_back(2,i,iset),qin_back(2,i,iset),
825 & wfrag_back(3,i,iset),qin_back(3,i,iset)
828 write(iout,*) "fragment, weights:"
830 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
831 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
832 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
839 call rescale_weights(t_bath)
841 print *,"chainbuild_cart"
843 print *,"After cartprint"
846 print *,"before ETOT"
847 write (iout,*) "usampl",usampl
848 call etotal(energy(0))
850 call enerprint(energy(0))
851 write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
852 print *,'icheckgrad=',icheckgrad
854 goto (10,20,30) icheckgrad
855 10 call check_ecartint
858 & "Checking the gradient of Cartesian coordinates disabled."
863 c---------------------------------------------------------------------------
870 c---------------------------------------------------------------------------
877 include 'COMMON.IOUNITS'
878 C Conformational Space Annealling programmed by Jooyoung Lee.
879 C This method works only with parallel machines!
883 write (iout,*) "CSA works on parallel machines only"
887 c---------------------------------------------------------------------------
888 subroutine exec_softreg
891 include 'COMMON.IOUNITS'
892 include 'COMMON.CONTROL'
893 double precision energy(0:max_ene),etot
894 double precision rms,frac,frac_nn,co
896 call etotal(energy(0))
897 call enerprint(energy(0))
898 if (.not.lsecondary) then
899 write(iout,*) 'Calling secondary structure recognition'
900 call secondary2(.true.)
902 write(iout,*) 'Using secondary structure supplied in pdb'
907 call etotal(energy(0))
909 call enerprint(energy(0))
911 call briefout(0,etot)
912 call secondary2(.true.)
913 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)