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,it,icall,iretcode,nfun
218 integer nharp,iharp(4,maxres/3)
220 logical fail,secondary_str /.true./
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)
246 if (nran_start.gt.0) then
248 & "Chains will be regenerated starting from residue",nran_start
250 call gen_rand_conf_mchain(nran_start,*10)
251 write (iout,*) "Conformation successfully generated",it
253 10 write (iout,*) "Problems with regenerating chains",it
256 write (iout,*) "Cartesian coords after chain rebuild"
259 write (iout,*) "Cartesian coords after chainbuild_ecart"
261 call int_from_cart1(.false.)
265 print *,"Processor",myrank," after chainbuild"
267 call etotal_long(energy_long(0))
268 write (iout,*) "Printing long range energy"
269 call enerprint(energy_long(0))
270 call etotal_short(energy_short(0))
271 write (iout,*) "Printing short range energy"
272 call enerprint(energy_short(0))
274 energy(i)=energy_long(i)+energy_short(i)
275 c write (iout,*) i,energy_long(i),energy_short(i),energy(i)
277 write (iout,*) "Printing long+short range energy"
278 call enerprint(energy(0))
280 c write(iout,*)"before etotal"
282 call etotal(energy(0))
283 c write(iout,*)"after etotal"
286 time_ene=MPI_Wtime()-time00
288 time_ene=tcpu()-time00
290 write (iout,*) "Time for energy evaluation",time_ene
291 c print *,"after etotal"
294 call enerprint(energy(0))
295 if (secondary_str) then
296 call hairpin(.true.,nharp,iharp)
297 c print *,'after hairpin'
298 call secondary2(.true.)
300 c print *,'after secondary'
303 if (indpdb.ne.0 .and. .not.dccart) then
305 call chainbuild_extconf
306 call etotal(energy(0))
307 write (iout,*) "After bond regularization"
308 call enerprint(energy(0))
312 write (iout,*) 'Calling OVERLAP_SC'
313 call overlap_sc(fail)
314 write (iout,*) "After overlap_sc"
315 c cartname=prefix(:ilen(prefix))//'.x'
317 c call cartoutx(0.0d0)
321 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
322 print *,'SC_move',nft_sc,etot
323 write(iout,*) 'SC_move',nft_sc,etot
327 print *, 'Calling MINIM_DC'
333 call minim_dc(etot,iretcode,nfun)
335 call geom_to_var(nvar,varia)
336 c print *,'Calling MINIMIZE.'
342 call minimize(etot,varia,iretcode,nfun)
345 print *,'LBFGS return code is',status,' eval ',nfun
347 print *,'SUMSL return code is',iretcode,' eval ',nfun
350 evals=nfun/(MPI_WTIME()-time1)
352 evals=nfun/(tcpu()-time1)
354 print *,'# eval/s',evals
355 print *,'refstr=',refstr
356 if (secondary_str) then
357 call hairpin(.false.,nharp,iharp)
358 c print *,'after hairpin'
359 call secondary2(.true.)
360 c print *,'after secondary'
362 call etotal(energy(0))
364 call enerprint(energy(0))
367 write (iout,'(a,a9)') 'LBFGS return code:',status
368 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
369 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
371 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
372 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
373 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
376 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
377 if (out_int) call briefout(0,etot)
379 cartname=prefix(:ilen(prefix))//'.x'
383 if (outpdb) call pdbout(etot,titel(:50),ipdb)
384 if (outmol2) call mol2out(etot,titel(:32))
387 c---------------------------------------------------------------------------
388 subroutine exec_regularize
394 include 'COMMON.SETUP'
395 include 'COMMON.TIME1'
396 include 'COMMON.INTERACT'
397 include 'COMMON.NAMES'
399 include 'COMMON.HEADER'
400 include 'COMMON.CONTROL'
401 c include 'COMMON.CONTACTS'
402 include 'COMMON.CHAIN'
404 include 'COMMON.IOUNITS'
405 include 'COMMON.FFIELD'
406 include 'COMMON.REMD'
408 include 'COMMON.SBRIDGE'
409 double precision energy(0:n_ene)
410 double precision etot,rms,frac,frac_nn,co
415 common /lbfgstat/ status,niter,nfun
421 call regularize(nct-nnt+1,etot,rms,cref(1,nnt),iretcode)
422 call etotal(energy(0))
423 energy(0)=energy(0)-energy(14)
425 call enerprint(energy(0))
427 call briefout(0,etot)
428 if (outpdb) call pdbout(etot,titel(:50),ipdb)
429 if (outmol2) call mol2out(etot,titel(:32))
430 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
432 write (iout,'(a,a9)') 'LBFGS return code:',status
434 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
438 c---------------------------------------------------------------------------
439 subroutine exec_thread
445 include "COMMON.SETUP"
449 c---------------------------------------------------------------------------
453 character*10 nodeinfo
455 double precision varia(maxvar)
459 include "COMMON.SETUP"
460 include 'COMMON.CONTROL'
464 if (modecalc.eq.3) then
470 if (modecalc.eq.3) then
481 c---------------------------------------------------------------------------
482 subroutine exec_mult_eeval_or_minim
487 integer muster(mpi_status_size)
490 include 'COMMON.SETUP'
491 include 'COMMON.TIME1'
492 include 'COMMON.INTERACT'
493 include 'COMMON.NAMES'
495 include 'COMMON.HEADER'
496 include 'COMMON.CONTROL'
497 c include 'COMMON.CONTACTS'
498 include 'COMMON.CHAIN'
500 include 'COMMON.IOUNITS'
501 include 'COMMON.FFIELD'
502 include 'COMMON.REMD'
504 include 'COMMON.SBRIDGE'
505 double precision varia(maxvar)
506 integer i,j,iconf,ind(6)
507 integer n,it,man,nf_mcmf,nmin,imm,mm,nft
508 double precision energy(0:max_ene),ene,etot,ene0
509 double precision rms,frac,frac_nn,co
510 double precision time
520 open(intin,file=intinname,status='old')
521 write (istat,'(a5,20a12)')"# ",
522 & (wname(print_order(i)),i=1,nprint_ene)
524 write (istat,'(a5,20a12)')"# ",
525 & (ename(print_order(i)),i=1,nprint_ene),
526 & "ETOT total","RMSD","nat.contact","nnt.contact"
528 write (istat,'(a5,20a12)')"# ",
529 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
535 read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
536 call read_x(intin,*11)
538 c Broadcast the order to compute internal coordinates to the slaves.
540 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
542 call int_from_cart1(.false.)
544 read (intin,'(i5)',end=1100,err=1100) iconf
545 call read_angles(intin,*11)
546 call geom_to_var(nvar,varia)
549 write (iout,'(a,i7)') 'Conformation #',iconf
550 call etotal(energy(0))
551 call briefout(iconf,energy(0))
552 call enerprint(energy(0))
555 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
556 write (istat,'(i5,20(f12.3))') iconf,
557 & (energy(print_order(i)),i=1,nprint_ene),etot,
558 & rms,frac,frac_nn,co
561 write (istat,'(i5,16(f12.3))') iconf,
562 & (energy(print_order(i)),i=1,nprint_ene),etot
578 if (mm.lt.nodes) then
580 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
581 call read_x(intin,*11)
583 c Broadcast the order to compute internal coordinates to the slaves.
585 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
587 call int_from_cart1(.false.)
589 read (intin,'(i5)',end=11,err=11) iconf
590 call read_angles(intin,*11)
591 call geom_to_var(nvar,varia)
594 write (iout,'(a,i7)') 'Conformation #',iconf
604 call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
606 call mpi_send(varia,nvar,mpi_double_precision,mm,
607 * idreal,CG_COMM,ierr)
608 call mpi_send(ene0,1,mpi_double_precision,mm,
609 * idreal,CG_COMM,ierr)
610 c print *,'task ',n,' sent to worker ',mm,nvar
612 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
613 * CG_COMM,muster,ierr)
614 man=muster(mpi_source)
615 c print *,'receiving result from worker ',man,' (',iii1,iii,')'
616 call mpi_recv(varia,nvar,mpi_double_precision,
617 * man,idreal,CG_COMM,muster,ierr)
619 * mpi_double_precision,man,idreal,
620 * CG_COMM,muster,ierr)
621 call mpi_recv(ene0,1,
622 * mpi_double_precision,man,idreal,
623 * CG_COMM,muster,ierr)
624 c print *,'result received from worker ',man,' sending now'
626 call var_to_geom(nvar,varia)
628 call etotal(energy(0))
632 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
635 call enerprint(energy(0))
636 call briefout(it,etot)
637 c if (minim) call briefout(it,etot)
639 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
640 write (istat,'(i5,19(f12.3))') iconf,
641 & (energy(print_order(i)),i=1,nprint_ene),etot,
642 & rms,frac,frac_nn,co
644 write (istat,'(i5,15(f12.3))') iconf,
645 & (energy(print_order(i)),i=1,nprint_ene),etot
650 read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
651 call read_x(intin,*11)
653 c Broadcast the order to compute internal coordinates to the slaves.
655 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
657 call int_from_cart1(.false.)
659 read (intin,'(i5)',end=1101,err=1101) iconf
660 call read_angles(intin,*11)
661 call geom_to_var(nvar,varia)
672 call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
674 call mpi_send(varia,nvar,mpi_double_precision,man,
675 * idreal,CG_COMM,ierr)
676 call mpi_send(ene0,1,mpi_double_precision,man,
677 * idreal,CG_COMM,ierr)
678 nf_mcmf=nf_mcmf+ind(4)
684 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
685 * CG_COMM,muster,ierr)
686 man=muster(mpi_source)
687 call mpi_recv(varia,nvar,mpi_double_precision,
688 * man,idreal,CG_COMM,muster,ierr)
690 * mpi_double_precision,man,idreal,
691 * CG_COMM,muster,ierr)
692 call mpi_recv(ene0,1,
693 * mpi_double_precision,man,idreal,
694 * CG_COMM,muster,ierr)
696 call var_to_geom(nvar,varia)
698 call etotal(energy(0))
702 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
705 call enerprint(energy(0))
706 call briefout(it,etot)
708 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
709 write (istat,'(i5,19(f12.3))') iconf,
710 & (energy(print_order(i)),i=1,nprint_ene),etot,
711 & rms,frac,frac_nn,co
713 write (istat,'(i5,15(f12.3))') iconf,
714 & (energy(print_order(i)),i=1,nprint_ene),etot
726 call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
731 open(intin,file=intinname,status='old')
732 write (istat,'(a5,20a12)')"# ",
733 & (wname(print_order(i)),i=1,nprint_ene)
734 write (istat,'("# ",20(1pe12.4))')
735 & (weights(print_order(i)),i=1,nprint_ene)
737 write (istat,'(a5,20a12)')"# ",
738 & (ename(print_order(i)),i=1,nprint_ene),
739 & "ETOT total","RMSD","nat.contact","nnt.contact"
741 write (istat,'(a5,14a12)')"# ",
742 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
746 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
747 call read_x(intin,*11)
749 c Broadcast the order to compute internal coordinates to the slaves.
751 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
753 call int_from_cart1(.false.)
755 read (intin,'(i5)',end=11,err=11) iconf
756 call read_angles(intin,*11)
757 call geom_to_var(nvar,varia)
760 write (iout,'(a,i7)') 'Conformation #',iconf
761 if (minim) call minimize(etot,varia,iretcode,nfun)
762 call etotal(energy(0))
765 call enerprint(energy(0))
766 if (minim) call briefout(it,etot)
768 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
769 write (istat,'(i5,18(f12.3))') iconf,
770 & (energy(print_order(i)),i=1,nprint_ene),
771 & etot,rms,frac,frac_nn,co
774 write (istat,'(i5,14(f12.3))') iconf,
775 & (energy(print_order(i)),i=1,nprint_ene),etot
782 c---------------------------------------------------------------------------
783 subroutine exec_checkgrad
789 include 'COMMON.SETUP'
790 include 'COMMON.TIME1'
791 include 'COMMON.INTERACT'
792 include 'COMMON.NAMES'
794 include 'COMMON.HEADER'
795 include 'COMMON.CONTROL'
796 c include 'COMMON.CONTACTS'
797 include 'COMMON.CHAIN'
799 include 'COMMON.IOUNITS'
800 include 'COMMON.FFIELD'
801 include 'COMMON.REMD'
803 include 'COMMON.QRESTR'
804 include 'COMMON.SBRIDGE'
807 double precision energy(0:max_ene)
810 c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
811 c if (itype(i).ne.10)
812 c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
814 if (indpdb.eq.0) call chainbuild
817 c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
821 c if (itype(i).ne.10) then
823 c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
828 c dc(j,0)=ran_number(-0.2d0,0.2d0)
841 print *, "AFTER read fragments"
842 write (iout,*) "iset",iset
844 write(iout,*) "fragment, weights, q0:"
846 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
847 & ifrag_back(2,i,iset),
848 & wfrag_back(1,i,iset),qin_back(1,i,iset),
849 & wfrag_back(2,i,iset),qin_back(2,i,iset),
850 & wfrag_back(3,i,iset),qin_back(3,i,iset)
853 write(iout,*) "fragment, weights:"
855 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
856 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
857 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
864 call rescale_weights(t_bath)
866 print *,"chainbuild_cart"
868 print *,"After cartprint"
871 print *,"before ETOT"
872 write (iout,*) "usampl",usampl
873 call etotal(energy(0))
875 call enerprint(energy(0))
876 write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
877 print *,'icheckgrad=',icheckgrad
879 goto (10,20,30) icheckgrad
880 10 call check_ecartint
883 & "Checking the gradient of Cartesian coordinates disabled."
888 c---------------------------------------------------------------------------
895 c---------------------------------------------------------------------------
902 include 'COMMON.IOUNITS'
903 C Conformational Space Annealling programmed by Jooyoung Lee.
904 C This method works only with parallel machines!
908 write (iout,*) "CSA works on parallel machines only"
912 c---------------------------------------------------------------------------
913 subroutine exec_softreg
916 include 'COMMON.IOUNITS'
917 include 'COMMON.CONTROL'
918 double precision energy(0:max_ene),etot
919 double precision rms,frac,frac_nn,co
921 call etotal(energy(0))
922 call enerprint(energy(0))
923 if (.not.lsecondary) then
924 write(iout,*) 'Calling secondary structure recognition'
925 call secondary2(.true.)
927 write(iout,*) 'Using secondary structure supplied in pdb'
932 call etotal(energy(0))
934 call enerprint(energy(0))
936 call briefout(0,etot)
937 call secondary2(.true.)
938 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)