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 write (iout,*) "After readrtns"
66 if (me.eq.king .or. .not. out1file) then
68 & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
70 if (minim) write (iout,'(a)')
71 & 'Conformations will be energy-minimized.'
72 write (iout,'(80(1h*)/)')
76 if (modecalc.eq.-2) then
79 else if (modecalc.eq.-1) then
80 write(iout,*) "call check_sc_map next"
85 if (fg_rank.gt.0) then
86 C Fine-grain slaves just do energy and gradient components.
87 call ergastulum ! slave workhouse in Latin
90 if (indpdb.eq.0 .and. .not.read_cart) call chainbuild
91 if (indpdb.ne.0 .or. read_cart) then
96 if (modecalc.eq.0) then
97 write (iout,*) "Calling exec_eeval_or_minim"
98 call exec_eeval_or_minim
99 else if (modecalc.eq.1) then
101 else if (modecalc.eq.2) then
103 else if (modecalc.eq.3 .or. modecalc .eq.6) then
105 else if (modecalc.eq.4) then
106 call exec_mult_eeval_or_minim
107 else if (modecalc.eq.5) then
109 else if (ModeCalc.eq.7) then
111 else if (ModeCalc.eq.8) then
113 else if (modecalc.eq.11) then
115 else if (modecalc.eq.12) then
117 else if (modecalc.eq.14) then
121 write (iout,*) "Need a parallel version to run MREMD."
125 write (iout,'(a)') 'This calculation type is not supported',
131 if (fg_rank.eq.0) call finish_task
132 c call memmon_print_usage()
134 call print_detailed_timing
136 call MPI_Finalize(ierr)
139 call dajczas(tcpu(),hrtime,mintime,sectime)
140 stop '********** Program terminated normally.'
143 c--------------------------------------------------------------------------
150 include 'COMMON.SETUP'
151 include 'COMMON.CONTROL'
152 include 'COMMON.IOUNITS'
153 c if (me.eq.king .or. .not. out1file) then
154 c write (iout,*) "Calling chainbuild"
158 c if (me.eq.king .or. .not. out1file) then
159 c write (iout,*) "Calling MD"
165 c---------------------------------------------------------------------------
167 subroutine exec_MREMD
173 include 'COMMON.SETUP'
174 include 'COMMON.CONTROL'
175 include 'COMMON.IOUNITS'
176 include 'COMMON.REMD'
178 if (me.eq.king .or. .not. out1file)
179 & write (iout,*) "Calling chainbuild"
181 if (me.eq.king .or. .not. out1file)
182 & write (iout,*) "Calling REMD"
194 c---------------------------------------------------------------------------
195 subroutine exec_eeval_or_minim
201 include 'COMMON.SETUP'
202 include 'COMMON.TIME1'
203 include 'COMMON.INTERACT'
204 include 'COMMON.NAMES'
206 include 'COMMON.HEADER'
207 include 'COMMON.CONTROL'
208 c include 'COMMON.CONTACTS'
209 include 'COMMON.CHAIN'
211 include 'COMMON.IOUNITS'
212 include 'COMMON.FFIELD'
213 include 'COMMON.REMD'
215 include 'COMMON.SBRIDGE'
216 integer i,icall,iretcode,nfun
218 integer nharp,iharp(4,maxres/3)
221 double precision energy(0:n_ene),etot,etota
222 double precision energy_long(0:n_ene),energy_short(0:n_ene)
223 double precision rms,frac,frac_nn,co
224 double precision varia(maxvar)
225 double precision time00,time1,time_ene,evals
229 common /lbfgstat/ status,niter,nfun
232 if (indpdb.eq.0 .and. .not.read_cart) call chainbuild
233 if (indpdb.ne.0 .or. read_cart) then
243 write (iout,*) "Energy evaluation/minimization"
245 c print *,'dc',dc(1,0),dc(2,0),dc(3,0)
247 print *,"Processor",myrank," after chainbuild"
249 call etotal_long(energy_long(0))
250 write (iout,*) "Printing long range energy"
251 call enerprint(energy_long(0))
252 call etotal_short(energy_short(0))
253 write (iout,*) "Printing short range energy"
254 call enerprint(energy_short(0))
256 energy(i)=energy_long(i)+energy_short(i)
257 c write (iout,*) i,energy_long(i),energy_short(i),energy(i)
259 write (iout,*) "Printing long+short range energy"
260 call enerprint(energy(0))
262 c write(iout,*)"before etotal"
264 call etotal(energy(0))
265 c write(iout,*)"after etotal"
268 time_ene=MPI_Wtime()-time00
270 time_ene=tcpu()-time00
272 write (iout,*) "Time for energy evaluation",time_ene
273 c print *,"after etotal"
276 call enerprint(energy(0))
277 call hairpin(.true.,nharp,iharp)
278 c print *,'after hairpin'
279 call secondary2(.true.)
280 c print *,'after secondary'
283 if (indpdb.ne.0 .and. .not.dccart) then
285 call chainbuild_extconf
286 call etotal(energy(0))
287 write (iout,*) "After bond regularization"
288 call enerprint(energy(0))
292 write (iout,*) 'Calling OVERLAP_SC'
293 call overlap_sc(fail)
294 write (iout,*) "After overlap_sc"
298 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
299 print *,'SC_move',nft_sc,etot
300 write(iout,*) 'SC_move',nft_sc,etot
304 print *, 'Calling MINIM_DC'
310 call minim_dc(etot,iretcode,nfun)
312 call geom_to_var(nvar,varia)
313 c print *,'Calling MINIMIZE.'
319 call minimize(etot,varia,iretcode,nfun)
322 print *,'LBFGS return code is',status,' eval ',nfun
324 print *,'SUMSL return code is',iretcode,' eval ',nfun
327 evals=nfun/(MPI_WTIME()-time1)
329 evals=nfun/(tcpu()-time1)
331 print *,'# eval/s',evals
332 print *,'refstr=',refstr
333 call hairpin(.false.,nharp,iharp)
334 print *,'after hairpin'
335 call secondary2(.true.)
336 print *,'after secondary'
337 call etotal(energy(0))
339 call enerprint(energy(0))
342 write (iout,'(a,a9)') 'LBFGS return code:',status
343 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
344 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
346 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
347 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
348 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
351 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
352 if (out_int) call briefout(0,etot)
354 cartname=prefix(:ilen(prefix))//'.x'
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.)