X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Funres.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Funres.F;h=d5031f0b7cf416f2e9c28357c0f887ae7194898c;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/unres.F b/source/unres/src_MD-M-newcorr/unres.F new file mode 100644 index 0000000..d5031f0 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/unres.F @@ -0,0 +1,771 @@ +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C C +C U N R E S C +C C +C Program to carry out conformational search of proteins in an united-residue C +C approximation. C +C C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + + +#ifdef MPI + include 'mpif.h' + include 'COMMON.SETUP' +#endif + include 'COMMON.TIME1' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.GEO' + include 'COMMON.HEADER' + include 'COMMON.CONTROL' + include 'COMMON.CONTACTS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + include 'COMMON.MD' + include 'COMMON.SBRIDGE' + double precision hrtime,mintime,sectime + character*64 text_mode_calc(-2:14) /'test', + & 'SC rotamer distribution', + & 'Energy evaluation or minimization', + & 'Regularization of PDB structure', + & 'Threading of a sequence on PDB structures', + & 'Monte Carlo (with minimization) ', + & 'Energy minimization of multiple conformations', + & 'Checking energy gradient', + & 'Entropic sampling Monte Carlo (with minimization)', + & 'Energy map', + & 'CSA calculations', + & 'Not used 9', + & 'Not used 10', + & 'Soft regularization of PDB structure', + & 'Mesoscopic molecular dynamics (MD) ', + & 'Not used 13', + & 'Replica exchange molecular dynamics (REMD)'/ + external ilen + +c call memmon_print_usage() + + call init_task + if (me.eq.king) + & write(iout,*)'### LAST MODIFIED 4/25/08 7:29PM by adam' + if (me.eq.king) call cinfo +C Read force field parameters and job setup data + call readrtns + call flush(iout) +C + if (me.eq.king .or. .not. out1file) then + write (iout,'(2a/)') + & text_mode_calc(modecalc)(:ilen(text_mode_calc(modecalc))), + & ' calculation.' + if (minim) write (iout,'(a)') + & 'Conformations will be energy-minimized.' + write (iout,'(80(1h*)/)') + endif + call flush(iout) +C + if (modecalc.eq.-2) then + call test + stop + else if (modecalc.eq.-1) then + write(iout,*) "call check_sc_map next" + call check_bond + stop + endif +#ifdef MPI + if (fg_rank.gt.0) then +C Fine-grain slaves just do energy and gradient components. + call ergastulum ! slave workhouse in Latin + else +#endif + if (modecalc.eq.0) then + call exec_eeval_or_minim + else if (modecalc.eq.1) then + call exec_regularize + else if (modecalc.eq.2) then + call exec_thread + else if (modecalc.eq.3 .or. modecalc .eq.6) then + call exec_MC + else if (modecalc.eq.4) then + call exec_mult_eeval_or_minim + else if (modecalc.eq.5) then + call exec_checkgrad + call exec_checkgrad + else if (ModeCalc.eq.7) then + call exec_map + else if (ModeCalc.eq.8) then + call exec_CSA + else if (modecalc.eq.11) then + call exec_softreg + else if (modecalc.eq.12) then + call exec_MD + else if (modecalc.eq.14) then + call exec_MREMD + else + write (iout,'(a)') 'This calculation type is not supported', + & ModeCalc + endif +#ifdef MPI + endif +C Finish task. + if (fg_rank.eq.0) call finish_task +c call memmon_print_usage() +#ifdef TIMING + call print_detailed_timing +#endif + call MPI_Finalize(ierr) + stop 'Bye Bye...' +#else + call dajczas(tcpu(),hrtime,mintime,sectime) + stop '********** Program terminated normally.' +#endif + end +c-------------------------------------------------------------------------- + subroutine exec_MD + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' + if (me.eq.king .or. .not. out1file) + & write (iout,*) "Calling chainbuild" + call chainbuild + call MD + return + end +c--------------------------------------------------------------------------- + subroutine exec_MREMD + include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' + include 'COMMON.REMD' + if (me.eq.king .or. .not. out1file) + & write (iout,*) "Calling chainbuild" + call chainbuild + if (me.eq.king .or. .not. out1file) + & write (iout,*) "Calling REMD" + if (remd_mlist) then + call MREMD + else + do i=1,nrep + remd_m(i)=1 + enddo + call MREMD + endif + return + end +c--------------------------------------------------------------------------- + subroutine exec_eeval_or_minim + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.TIME1' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.GEO' + include 'COMMON.HEADER' + include 'COMMON.CONTROL' + include 'COMMON.CONTACTS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + include 'COMMON.MD' + include 'COMMON.SBRIDGE' + common /srutu/ icall + double precision energy(0:n_ene) + double precision energy_long(0:n_ene),energy_short(0:n_ene) + double precision varia(maxvar) + if (indpdb.eq.0) call chainbuild + time00=MPI_Wtime() + call chainbuild_cart + if (split_ene) then + print *,"Processor",myrank," after chainbuild" + icall=1 + call etotal_long(energy_long(0)) + write (iout,*) "Printing long range energy" + call enerprint(energy_long(0)) + call etotal_short(energy_short(0)) + write (iout,*) "Printing short range energy" + call enerprint(energy_short(0)) + do i=0,n_ene + energy(i)=energy_long(i)+energy_short(i) + write (iout,*) i,energy_long(i),energy_short(i),energy(i) + enddo + write (iout,*) "Printing long+short range energy" + call enerprint(energy(0)) + endif + call etotal(energy(0)) + time_ene=MPI_Wtime()-time00 + write (iout,*) "Time for energy evaluation",time_ene + print *,"after etotal" + etota = energy(0) + etot =etota + call enerprint(energy(0)) + call hairpin(.true.,nharp,iharp) + call secondary2(.true.) + if (minim) then +crc overlap test + if (overlapsc) then + print *, 'Calling OVERLAP_SC' + call overlap_sc(fail) + endif + + if (searchsc) then + call sc_move(2,nres-1,10,1d10,nft_sc,etot) + print *,'SC_move',nft_sc,etot + write(iout,*) 'SC_move',nft_sc,etot + endif + + if (dccart) then + print *, 'Calling MINIM_DC' + time1=MPI_WTIME() + call minim_dc(etot,iretcode,nfun) + else + if (indpdb.ne.0) then + call bond_regular + call chainbuild + endif + call geom_to_var(nvar,varia) + print *,'Calling MINIMIZE.' + time1=MPI_WTIME() + call minimize(etot,varia,iretcode,nfun) + endif + print *,'SUMSL return code is',iretcode,' eval ',nfun + evals=nfun/(MPI_WTIME()-time1) + print *,'# eval/s',evals + print *,'refstr=',refstr + call hairpin(.true.,nharp,iharp) + call secondary2(.true.) + call etotal(energy(0)) + etot = energy(0) + call enerprint(energy(0)) + + call intout + call briefout(0,etot) + if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (iout,'(a,i3)') 'SUMSL return code:',iretcode + write (iout,'(a,i20)') '# of energy evaluations:',nfun+1 + write (iout,'(a,f16.3)')'# of energy evaluations/sec:',evals + else + print *,'refstr=',refstr + if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + call briefout(0,etot) + endif + if (outpdb) call pdbout(etot,titel(:50),ipdb) + if (outmol2) call mol2out(etot,titel) + return + end +c--------------------------------------------------------------------------- + subroutine exec_regularize + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.TIME1' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.GEO' + include 'COMMON.HEADER' + include 'COMMON.CONTROL' + include 'COMMON.CONTACTS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + include 'COMMON.MD' + include 'COMMON.SBRIDGE' + double precision energy(0:n_ene) + + call gen_dist_constr + call sc_conf + call intout + call regularize(nct-nnt+1,etot,rms,cref(1,nnt,1),iretcode) + call etotal(energy(0)) + energy(0)=energy(0)-energy(14) + etot=energy(0) + call enerprint(energy(0)) + call intout + call briefout(0,etot) + if (outpdb) call pdbout(etot,titel,ipdb) + if (outmol2) call mol2out(etot,titel(:32)) + if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (iout,'(a,i3)') 'SUMSL return code:',iretcode + return + end +c--------------------------------------------------------------------------- + subroutine exec_thread + include 'DIMENSIONS' +#ifdef MP + include "mpif.h" +#endif + include "COMMON.SETUP" + call thread_seq + return + end +c--------------------------------------------------------------------------- + subroutine exec_MC + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + character*10 nodeinfo + double precision varia(maxvar) +#ifdef MPI + include "mpif.h" +#endif + include "COMMON.SETUP" + include 'COMMON.CONTROL' + call mcm_setup + if (minim) then +#ifdef MPI + if (modecalc.eq.3) then + call do_mcm(ipar) + else + call entmcm + endif +#else + if (modecalc.eq.3) then + call do_mcm(ipar) + else + call entmcm + endif +#endif + else + call monte_carlo + endif + return + end +c--------------------------------------------------------------------------- + subroutine exec_mult_eeval_or_minim + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' + dimension muster(mpi_status_size) +#endif + include 'COMMON.SETUP' + include 'COMMON.TIME1' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.GEO' + include 'COMMON.HEADER' + include 'COMMON.CONTROL' + include 'COMMON.CONTACTS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + include 'COMMON.MD' + include 'COMMON.SBRIDGE' + double precision varia(maxvar) + dimension ind(6) + double precision energy(0:max_ene) + logical eof + eof=.false. +#ifdef MPI + if(me.ne.king) then + call minim_mcmf + return + endif + + close (intin) + open(intin,file=intinname,status='old') + write (istat,'(a5,20a12)')"# ", + & (wname(print_order(i)),i=1,nprint_ene) + if (refstr) then + write (istat,'(a5,20a12)')"# ", + & (ename(print_order(i)),i=1,nprint_ene), + & "ETOT total","RMSD","nat.contact","nnt.contact" + else + write (istat,'(a5,20a12)')"# ", + & (ename(print_order(i)),i=1,nprint_ene),"ETOT total" + endif + + if (.not.minim) then + do while (.not. eof) + if (read_cart) then + read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene + call read_x(intin,*11) +#ifdef MPI +c Broadcast the order to compute internal coordinates to the slaves. + if (nfgtasks.gt.1) + & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) +#endif + call int_from_cart1(.false.) + else + read (intin,'(i5)',end=1100,err=1100) iconf + call read_angles(intin,*11) + call geom_to_var(nvar,varia) + call chainbuild + endif + write (iout,'(a,i7)') 'Conformation #',iconf + call etotal(energy(0)) + call briefout(iconf,energy(0)) + call enerprint(energy(0)) + etot=energy(0) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (istat,'(i5,20(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot, + & rms,frac,frac_nn,co +cjlee end + else + write (istat,'(i5,16(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot + endif + enddo +1100 continue + goto 1101 + endif + + mm=0 + imm=0 + nft=0 + ene0=0.0d0 + n=0 + iconf=0 +c do n=1,nzsc + do while (.not. eof) + mm=mm+1 + if (mm.lt.nodes) then + if (read_cart) then + read (intin,'(e15.10,e15.5)',end=11,err=11) time,ene + call read_x(intin,*11) +#ifdef MPI +c Broadcast the order to compute internal coordinates to the slaves. + if (nfgtasks.gt.1) + & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) +#endif + call int_from_cart1(.false.) + else + read (intin,'(i5)',end=11,err=11) iconf + call read_angles(intin,*11) + call geom_to_var(nvar,varia) + call chainbuild + endif + write (iout,'(a,i7)') 'Conformation #',iconf + n=n+1 + imm=imm+1 + ind(1)=1 + ind(2)=n + ind(3)=0 + ind(4)=0 + ind(5)=0 + ind(6)=0 + ene0=0.0d0 + call mpi_send(ind,6,mpi_integer,mm,idint,CG_COMM, + * ierr) + call mpi_send(varia,nvar,mpi_double_precision,mm, + * idreal,CG_COMM,ierr) + call mpi_send(ene0,1,mpi_double_precision,mm, + * idreal,CG_COMM,ierr) +c print *,'task ',n,' sent to worker ',mm,nvar + else + call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint, + * CG_COMM,muster,ierr) + man=muster(mpi_source) +c print *,'receiving result from worker ',man,' (',iii1,iii,')' + call mpi_recv(varia,nvar,mpi_double_precision, + * man,idreal,CG_COMM,muster,ierr) + call mpi_recv(ene,1, + * mpi_double_precision,man,idreal, + * CG_COMM,muster,ierr) + call mpi_recv(ene0,1, + * mpi_double_precision,man,idreal, + * CG_COMM,muster,ierr) +c print *,'result received from worker ',man,' sending now' + + call var_to_geom(nvar,varia) + call chainbuild + call etotal(energy(0)) + iconf=ind(2) + write (iout,*) + write (iout,*) + write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5) + + etot=energy(0) + call enerprint(energy(0)) + call briefout(it,etot) +c if (minim) call briefout(it,etot) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (istat,'(i5,19(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot, + & rms,frac,frac_nn,co + else + write (istat,'(i5,15(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot + endif + + imm=imm-1 + if (read_cart) then + read (intin,'(e15.10,e15.5)',end=1101,err=1101) time,ene + call read_x(intin,*11) +#ifdef MPI +c Broadcast the order to compute internal coordinates to the slaves. + if (nfgtasks.gt.1) + & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) +#endif + call int_from_cart1(.false.) + else + read (intin,'(i5)',end=1101,err=1101) iconf + call read_angles(intin,*11) + call geom_to_var(nvar,varia) + call chainbuild + endif + n=n+1 + imm=imm+1 + ind(1)=1 + ind(2)=n + ind(3)=0 + ind(4)=0 + ind(5)=0 + ind(6)=0 + call mpi_send(ind,6,mpi_integer,man,idint,CG_COMM, + * ierr) + call mpi_send(varia,nvar,mpi_double_precision,man, + * idreal,CG_COMM,ierr) + call mpi_send(ene0,1,mpi_double_precision,man, + * idreal,CG_COMM,ierr) + nf_mcmf=nf_mcmf+ind(4) + nmin=nmin+1 + endif + enddo +11 continue + do j=1,imm + call mpi_recv(ind,6,mpi_integer,mpi_any_source,idint, + * CG_COMM,muster,ierr) + man=muster(mpi_source) + call mpi_recv(varia,nvar,mpi_double_precision, + * man,idreal,CG_COMM,muster,ierr) + call mpi_recv(ene,1, + * mpi_double_precision,man,idreal, + * CG_COMM,muster,ierr) + call mpi_recv(ene0,1, + * mpi_double_precision,man,idreal, + * CG_COMM,muster,ierr) + + call var_to_geom(nvar,varia) + call chainbuild + call etotal(energy(0)) + iconf=ind(2) + write (iout,*) + write (iout,*) + write (iout,'(a,2i7)') 'Conformation #',iconf,ind(5) + + etot=energy(0) + call enerprint(energy(0)) + call briefout(it,etot) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (istat,'(i5,19(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot, + & rms,frac,frac_nn,co + else + write (istat,'(i5,15(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot + endif + nmin=nmin+1 + enddo +1101 continue + do i=1, nodes-1 + ind(1)=0 + ind(2)=0 + ind(3)=0 + ind(4)=0 + ind(5)=0 + ind(6)=0 + call mpi_send(ind,6,mpi_integer,i,idint,CG_COMM, + * ierr) + enddo +#else + close (intin) + open(intin,file=intinname,status='old') + write (istat,'(a5,20a12)')"# ", + & (wname(print_order(i)),i=1,nprint_ene) + write (istat,'("# ",20(1pe12.4))') + & (weights(print_order(i)),i=1,nprint_ene) + if (refstr) then + write (istat,'(a5,20a12)')"# ", + & (ename(print_order(i)),i=1,nprint_ene), + & "ETOT total","RMSD","nat.contact","nnt.contact" + else + write (istat,'(a5,14a12)')"# ", + & (ename(print_order(i)),i=1,nprint_ene),"ETOT total" + endif + do while (.not. eof) + if (read_cart) then + read (intin,'(e15.10,e15.5)',end=1100,err=1100) time,ene + call read_x(intin,*11) +#ifdef MPI +c Broadcast the order to compute internal coordinates to the slaves. + if (nfgtasks.gt.1) + & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) +#endif + call int_from_cart1(.false.) + else + read (intin,'(i5)',end=1100,err=1100) iconf + call read_angles(intin,*11) + call geom_to_var(nvar,varia) + call chainbuild + endif + write (iout,'(a,i7)') 'Conformation #',iconf + if (minim) call minimize(etot,varia,iretcode,nfun) + call etotal(energy(0)) + + etot=energy(0) + call enerprint(energy(0)) + if (minim) call briefout(it,etot) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + write (istat,'(i5,18(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene), + & etot,rms,frac,frac_nn,co +cjlee end + else + write (istat,'(i5,14(f12.3))') iconf, + & (energy(print_order(i)),i=1,nprint_ene),etot + endif + enddo + 11 continue +#endif + return + end +c--------------------------------------------------------------------------- + subroutine exec_checkgrad + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.TIME1' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.GEO' + include 'COMMON.HEADER' + include 'COMMON.CONTROL' + include 'COMMON.CONTACTS' + include 'COMMON.CHAIN' + include 'COMMON.VAR' + include 'COMMON.IOUNITS' + include 'COMMON.FFIELD' + include 'COMMON.REMD' + include 'COMMON.MD' + include 'COMMON.SBRIDGE' + common /srutu/ icall + double precision energy(0:max_ene) +c do i=2,nres +c vbld(i)=vbld(i)+ran_number(-0.1d0,0.1d0) +c if (itype(i).ne.10) +c & vbld(i+nres)=vbld(i+nres)+ran_number(-0.001d0,0.001d0) +c enddo + if (indpdb.eq.0) call chainbuild +c do i=0,nres +c do j=1,3 +c dc(j,i)=dc(j,i)+ran_number(-0.2d0,0.2d0) +c enddo +c enddo +c do i=1,nres-1 +c if (itype(i).ne.10) then +c do j=1,3 +c dc(j,i+nres)=dc(j,i+nres)+ran_number(-0.2d0,0.2d0) +c enddo +c endif +c enddo +c do j=1,3 +c dc(j,0)=ran_number(-0.2d0,0.2d0) +c enddo + usampl=.true. + totT=1.d0 + eq_time=0.0d0 + call read_fragments + call chainbuild_cart + call cartprint + call intout + icall=1 + call etotal(energy(0)) + etot = energy(0) + call enerprint(energy(0)) + write (iout,*) "Uconst",Uconst," Uconst_back",uconst_back + print *,'icheckgrad=',icheckgrad + goto (10,20,30) icheckgrad + 10 print *,'kurwa0' + call check_ecartint + print *,'kurwa' + call check_ecartint + return + 20 print *,'ja pierdole' + call check_cartgrad + call check_cartgrad + return + 30 call check_eint + return + end +c--------------------------------------------------------------------------- + subroutine exec_map +C Energy maps + call map_read + call map + return + end +c--------------------------------------------------------------------------- + subroutine exec_CSA +#ifdef MPI + include "mpif.h" +#endif + include 'DIMENSIONS' + include 'COMMON.IOUNITS' +C Conformational Space Annealling programmed by Jooyoung Lee. +C This method works only with parallel machines! +#ifdef MPI + call together +#else + write (iout,*) "CSA works on parallel machines only" +#endif + return + end +c--------------------------------------------------------------------------- + subroutine exec_softreg + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CONTROL' + double precision energy(0:max_ene) + call chainbuild + call etotal(energy(0)) + call enerprint(energy(0)) + if (.not.lsecondary) then + write(iout,*) 'Calling secondary structure recognition' + call secondary2(debug) + else + write(iout,*) 'Using secondary structure supplied in pdb' + endif + + call softreg + + call etotal(energy(0)) + etot=energy(0) + call enerprint(energy(0)) + call intout + call briefout(0,etot) + call secondary2(.true.) + if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.) + return + end