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 write (iout,*) "After readrtns"
62 if (me.eq.king .or. .not. out1file) then
64 & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))),
66 if (minim) write (iout,'(a)')
67 & 'Conformations will be energy-minimized.'
68 write (iout,'(80(1h*)/)')
72 if (modecalc.eq.-2) then
75 else if (modecalc.eq.-1) then
76 write(iout,*) "call check_sc_map next"
81 if (fg_rank.gt.0) then
82 C Fine-grain slaves just do energy and gradient components.
83 call ergastulum ! slave workhouse in Latin
86 if (modecalc.eq.0) then
87 write (iout,*) "Calling exec_eeval_or_minim"
89 call exec_eeval_or_minim
90 else if (modecalc.eq.1) then
92 else if (modecalc.eq.2) then
94 else if (modecalc.eq.3 .or. modecalc .eq.6) then
96 else if (modecalc.eq.4) then
97 call exec_mult_eeval_or_minim
98 else if (modecalc.eq.5) then
100 else if (ModeCalc.eq.7) then
102 else if (ModeCalc.eq.8) then
104 else if (modecalc.eq.11) then
106 else if (modecalc.eq.12) then
108 else if (modecalc.eq.14) then
112 write (iout,*) "Need a parallel version to run MREMD."
116 write (iout,'(a)') 'This calculation type is not supported',
122 if (fg_rank.eq.0) call finish_task
123 c call memmon_print_usage()
125 call print_detailed_timing
127 call MPI_Finalize(ierr)
130 call dajczas(tcpu(),hrtime,mintime,sectime)
131 stop '********** Program terminated normally.'
134 c--------------------------------------------------------------------------
140 include 'COMMON.SETUP'
141 include 'COMMON.CONTROL'
142 include 'COMMON.IOUNITS'
143 c if (me.eq.king .or. .not. out1file) then
144 c write (iout,*) "Calling chainbuild"
148 c if (me.eq.king .or. .not. out1file) then
149 c write (iout,*) "Calling MD"
155 c---------------------------------------------------------------------------
157 subroutine exec_MREMD
162 include 'COMMON.SETUP'
163 include 'COMMON.CONTROL'
164 include 'COMMON.IOUNITS'
165 include 'COMMON.REMD'
166 if (me.eq.king .or. .not. out1file)
167 & write (iout,*) "Calling chainbuild"
169 if (me.eq.king .or. .not. out1file)
170 & write (iout,*) "Calling REMD"
182 c---------------------------------------------------------------------------
183 subroutine exec_eeval_or_minim
184 implicit real*8 (a-h,o-z)
189 include 'COMMON.SETUP'
190 include 'COMMON.TIME1'
191 include 'COMMON.INTERACT'
192 include 'COMMON.NAMES'
194 include 'COMMON.HEADER'
195 include 'COMMON.CONTROL'
196 include 'COMMON.CONTACTS'
197 include 'COMMON.CHAIN'
199 include 'COMMON.IOUNITS'
200 include 'COMMON.FFIELD'
201 include 'COMMON.REMD'
203 include 'COMMON.SBRIDGE'
205 double precision energy(0:n_ene)
206 double precision energy_long(0:n_ene),energy_short(0:n_ene)
207 double precision varia(maxvar)
208 if (indpdb.eq.0) call chainbuild
209 if (indpdb.ne.0) then
219 write (iout,*) "Energy evaluation/minimization"
221 c print *,'dc',dc(1,0),dc(2,0),dc(3,0)
223 print *,"Processor",myrank," after chainbuild"
225 call etotal_long(energy_long(0))
226 write (iout,*) "Printing long range energy"
227 call enerprint(energy_long(0))
228 call etotal_short(energy_short(0))
229 write (iout,*) "Printing short range energy"
230 call enerprint(energy_short(0))
232 energy(i)=energy_long(i)+energy_short(i)
233 c write (iout,*) i,energy_long(i),energy_short(i),energy(i)
235 write (iout,*) "Printing long+short range energy"
236 call enerprint(energy(0))
238 write(iout,*)"before etotal"
240 call etotal(energy(0))
241 write(iout,*)"after etotal"
244 time_ene=MPI_Wtime()-time00
246 time_ene=tcpu()-time00
248 write (iout,*) "Time for energy evaluation",time_ene
249 print *,"after etotal"
252 call enerprint(energy(0))
253 call hairpin(.true.,nharp,iharp)
254 print *,'after hairpin'
255 call secondary2(.true.)
256 print *,'after secondary'
260 print *, 'Calling OVERLAP_SC'
261 call overlap_sc(fail)
265 call sc_move(2,nres-1,10,1d10,nft_sc,etot)
266 print *,'SC_move',nft_sc,etot
267 write(iout,*) 'SC_move',nft_sc,etot
271 print *, 'Calling MINIM_DC'
277 call minim_dc(etot,iretcode,nfun)
279 if (indpdb.ne.0) then
283 call geom_to_var(nvar,varia)
284 print *,'Calling MINIMIZE.'
290 call minimize(etot,varia,iretcode,nfun)
292 print *,'SUMSL return code is',iretcode,' eval ',nfun
294 evals=nfun/(MPI_WTIME()-time1)
296 evals=nfun/(tcpu()-time1)
298 print *,'# eval/s',evals
299 print *,'refstr=',refstr
300 call hairpin(.false.,nharp,iharp)
301 print *,'after hairpin'
302 call secondary2(.true.)
303 print *,'after secondary'
304 call etotal(energy(0))
306 call enerprint(energy(0))
309 call briefout(0,etot)
310 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
311 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
312 write (iout,'(a,i20)') '# of energy evaluations:',nfun+1
313 write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals
315 print *,'refstr=',refstr
316 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
317 call briefout(0,etot)
319 if (outpdb) call pdbout(etot,titel(:50),ipdb)
320 if (outmol2) call mol2out(etot,titel(:32))
323 c---------------------------------------------------------------------------
324 subroutine exec_regularize
325 implicit real*8 (a-h,o-z)
330 include 'COMMON.SETUP'
331 include 'COMMON.TIME1'
332 include 'COMMON.INTERACT'
333 include 'COMMON.NAMES'
335 include 'COMMON.HEADER'
336 include 'COMMON.CONTROL'
337 include 'COMMON.CONTACTS'
338 include 'COMMON.CHAIN'
340 include 'COMMON.IOUNITS'
341 include 'COMMON.FFIELD'
342 include 'COMMON.REMD'
344 include 'COMMON.SBRIDGE'
345 double precision energy(0:n_ene)
350 call regularize(nct-nnt+1,etot,rms,cref(1,nnt,1),iretcode)
351 call etotal(energy(0))
352 energy(0)=energy(0)-energy(14)
354 call enerprint(energy(0))
356 call briefout(0,etot)
357 if (outpdb) call pdbout(etot,titel(:50),ipdb)
358 if (outmol2) call mol2out(etot,titel(:32))
359 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
360 write (iout,'(a,i3)') 'SUMSL return code:',iretcode
363 c---------------------------------------------------------------------------
364 subroutine exec_thread
369 include "COMMON.SETUP"
373 c---------------------------------------------------------------------------
375 implicit real*8 (a-h,o-z)
377 character*10 nodeinfo
378 double precision varia(maxvar)
382 include "COMMON.SETUP"
383 include 'COMMON.CONTROL'
387 if (modecalc.eq.3) then
393 if (modecalc.eq.3) then
404 c---------------------------------------------------------------------------
405 subroutine exec_mult_eeval_or_minim
406 implicit real*8 (a-h,o-z)
410 dimension muster(mpi_status_size)
412 include 'COMMON.SETUP'
413 include 'COMMON.TIME1'
414 include 'COMMON.INTERACT'
415 include 'COMMON.NAMES'
417 include 'COMMON.HEADER'
418 include 'COMMON.CONTROL'
419 include 'COMMON.CONTACTS'
420 include 'COMMON.CHAIN'
422 include 'COMMON.IOUNITS'
423 include 'COMMON.FFIELD'
424 include 'COMMON.REMD'
426 include 'COMMON.SBRIDGE'
427 double precision varia(maxvar)
429 double precision energy(0:max_ene)
439 open(intin,file=intinname,status='old')
440 write (istat,'(a5,20a12)')"# ",
441 & (wname(print_order(i)),i=1,nprint_ene)
443 write (istat,'(a5,20a12)')"# ",
444 & (ename(print_order(i)),i=1,nprint_ene),
445 & "ETOT total","RMSD","nat.contact","nnt.contact"
447 write (istat,'(a5,20a12)')"# ",
448 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
454 read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene
455 call read_x(intin,*11)
457 c Broadcast the order to compute internal coordinates to the slaves.
459 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
461 call int_from_cart1(.false.)
463 read (intin,'(i5)',end=1100,err=1100) iconf
464 call read_angles(intin,*11)
465 call geom_to_var(nvar,varia)
468 write (iout,'(a,i7)') 'Conformation #',iconf
469 call etotal(energy(0))
470 call briefout(iconf,energy(0))
471 call enerprint(energy(0))
474 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
475 write (istat,'(i5,20(f12.3))') iconf,
476 & (energy(print_order(i)),i=1,nprint_ene),etot,
477 & rms,frac,frac_nn,co
480 write (istat,'(i5,16(f12.3))') iconf,
481 & (energy(print_order(i)),i=1,nprint_ene),etot
497 if (mm.lt.nodes) then
499 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
500 call read_x(intin,*11)
502 c Broadcast the order to compute internal coordinates to the slaves.
504 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
506 call int_from_cart1(.false.)
508 read (intin,'(i5)',end=11,err=11) iconf
509 call read_angles(intin,*11)
510 call geom_to_var(nvar,varia)
513 write (iout,'(a,i7)') 'Conformation #',iconf
523 call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM,
525 call mpi_send(varia,nvar,mpi_double_precision,mm,
526 * idreal,CG_COMM,ierr)
527 call mpi_send(ene0,1,mpi_double_precision,mm,
528 * idreal,CG_COMM,ierr)
529 c print *,'task ',n,' sent to worker ',mm,nvar
531 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
532 * CG_COMM,muster,ierr)
533 man=muster(mpi_source)
534 c print *,'receiving result from worker ',man,' (',iii1,iii,')'
535 call mpi_recv(varia,nvar,mpi_double_precision,
536 * man,idreal,CG_COMM,muster,ierr)
538 * mpi_double_precision,man,idreal,
539 * CG_COMM,muster,ierr)
540 call mpi_recv(ene0,1,
541 * mpi_double_precision,man,idreal,
542 * CG_COMM,muster,ierr)
543 c print *,'result received from worker ',man,' sending now'
545 call var_to_geom(nvar,varia)
547 call etotal(energy(0))
551 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
554 call enerprint(energy(0))
555 call briefout(it,etot)
556 c if (minim) call briefout(it,etot)
558 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
559 write (istat,'(i5,19(f12.3))') iconf,
560 & (energy(print_order(i)),i=1,nprint_ene),etot,
561 & rms,frac,frac_nn,co
563 write (istat,'(i5,15(f12.3))') iconf,
564 & (energy(print_order(i)),i=1,nprint_ene),etot
569 read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene
570 call read_x(intin,*11)
572 c Broadcast the order to compute internal coordinates to the slaves.
574 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
576 call int_from_cart1(.false.)
578 read (intin,'(i5)',end=1101,err=1101) iconf
579 call read_angles(intin,*11)
580 call geom_to_var(nvar,varia)
591 call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM,
593 call mpi_send(varia,nvar,mpi_double_precision,man,
594 * idreal,CG_COMM,ierr)
595 call mpi_send(ene0,1,mpi_double_precision,man,
596 * idreal,CG_COMM,ierr)
597 nf_mcmf=nf_mcmf+ind(4)
603 call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint,
604 * CG_COMM,muster,ierr)
605 man=muster(mpi_source)
606 call mpi_recv(varia,nvar,mpi_double_precision,
607 * man,idreal,CG_COMM,muster,ierr)
609 * mpi_double_precision,man,idreal,
610 * CG_COMM,muster,ierr)
611 call mpi_recv(ene0,1,
612 * mpi_double_precision,man,idreal,
613 * CG_COMM,muster,ierr)
615 call var_to_geom(nvar,varia)
617 call etotal(energy(0))
621 write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5)
624 call enerprint(energy(0))
625 call briefout(it,etot)
627 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
628 write (istat,'(i5,19(f12.3))') iconf,
629 & (energy(print_order(i)),i=1,nprint_ene),etot,
630 & rms,frac,frac_nn,co
632 write (istat,'(i5,15(f12.3))') iconf,
633 & (energy(print_order(i)),i=1,nprint_ene),etot
645 call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM,
650 open(intin,file=intinname,status='old')
651 write (istat,'(a5,20a12)')"# ",
652 & (wname(print_order(i)),i=1,nprint_ene)
653 write (istat,'("# ",20(1pe12.4))')
654 & (weights(print_order(i)),i=1,nprint_ene)
656 write (istat,'(a5,20a12)')"# ",
657 & (ename(print_order(i)),i=1,nprint_ene),
658 & "ETOT total","RMSD","nat.contact","nnt.contact"
660 write (istat,'(a5,14a12)')"# ",
661 & (ename(print_order(i)),i=1,nprint_ene),"ETOT total"
665 read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene
666 call read_x(intin,*11)
668 c Broadcast the order to compute internal coordinates to the slaves.
670 & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
672 call int_from_cart1(.false.)
674 read (intin,'(i5)',end=11,err=11) iconf
675 call read_angles(intin,*11)
676 call geom_to_var(nvar,varia)
679 write (iout,'(a,i7)') 'Conformation #',iconf
680 if (minim) call minimize(etot,varia,iretcode,nfun)
681 call etotal(energy(0))
684 call enerprint(energy(0))
685 if (minim) call briefout(it,etot)
687 call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
688 write (istat,'(i5,18(f12.3))') iconf,
689 & (energy(print_order(i)),i=1,nprint_ene),
690 & etot,rms,frac,frac_nn,co
693 write (istat,'(i5,14(f12.3))') iconf,
694 & (energy(print_order(i)),i=1,nprint_ene),etot
701 c---------------------------------------------------------------------------
702 subroutine exec_checkgrad
703 implicit real*8 (a-h,o-z)
708 include 'COMMON.SETUP'
709 include 'COMMON.TIME1'
710 include 'COMMON.INTERACT'
711 include 'COMMON.NAMES'
713 include 'COMMON.HEADER'
714 include 'COMMON.CONTROL'
715 include 'COMMON.CONTACTS'
716 include 'COMMON.CHAIN'
718 include 'COMMON.IOUNITS'
719 include 'COMMON.FFIELD'
720 include 'COMMON.REMD'
722 include 'COMMON.SBRIDGE'
724 double precision energy(0:max_ene)
727 c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0)
728 c if (itype(i).ne.10)
729 c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0)
731 if (indpdb.eq.0) call chainbuild
734 c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0)
738 c if (itype(i).ne.10) then
740 c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0)
745 c dc(j,0)=ran_number(-0.2d0,0.2d0)
756 print *, "AFTER read fragments"
757 write (iout,*) "iset",iset
759 write(iout,*) "fragment, weights, q0:"
761 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
762 & ifrag_back(2,i,iset),
763 & wfrag_back(1,i,iset),qin_back(1,i,iset),
764 & wfrag_back(2,i,iset),qin_back(2,i,iset),
765 & wfrag_back(3,i,iset),qin_back(3,i,iset)
768 write(iout,*) "fragment, weights:"
770 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
771 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
772 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
777 call rescale_weights(t_bath)
779 print *,"chainbuild_cart"
781 print *,"After cartprint"
784 print *,"before ETOT"
785 call etotal(energy(0))
787 call enerprint(energy(0))
788 write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back
789 print *,'icheckgrad=',icheckgrad
791 goto (10,20,30) icheckgrad
792 10 call check_ecartint
795 20 call check_cartgrad
800 c---------------------------------------------------------------------------
807 c---------------------------------------------------------------------------
813 include 'COMMON.IOUNITS'
814 C Conformational Space Annealling programmed by Jooyoung Lee.
815 C This method works only with parallel machines!
819 write (iout,*) "CSA works on parallel machines only"
823 c---------------------------------------------------------------------------
824 subroutine exec_softreg
826 include 'COMMON.IOUNITS'
827 include 'COMMON.CONTROL'
828 double precision energy(0:max_ene)
830 call etotal(energy(0))
831 call enerprint(energy(0))
832 if (.not.lsecondary) then
833 write(iout,*) 'Calling secondary structure recognition'
834 call secondary2(debug)
836 write(iout,*) 'Using secondary structure supplied in pdb'
841 call etotal(energy(0))
843 call enerprint(energy(0))
845 call briefout(0,etot)
846 call secondary2(.true.)
847 if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)