From 43dfdc33f5b338546258b5a882326f97a76b5bd4 Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 16 Feb 2021 15:11:53 +0100 Subject: [PATCH] Adam's unres update --- source/unres/src-HCD-5D/COMMON.CHAIN | 10 +- source/unres/src-HCD-5D/COMMON.CONTROL | 6 +- source/unres/src-HCD-5D/COMMON.INTERACT | 41 +- source/unres/src-HCD-5D/COMMON.MD | 4 +- source/unres/src-HCD-5D/COMMON.TIME1 | 10 +- source/unres/src-HCD-5D/COMMON.TORCNSTR | 5 +- source/unres/src-HCD-5D/DIMENSIONS | 15 +- source/unres/src-HCD-5D/MD_A-MTS.F | 13 +- source/unres/src-HCD-5D/cartprint.f | 2 +- source/unres/src-HCD-5D/chain_symmetry.F | 23 +- source/unres/src-HCD-5D/chainbuild.F | 37 +- source/unres/src-HCD-5D/checkder_p.F | 1 + source/unres/src-HCD-5D/contact.f | 6 +- source/unres/src-HCD-5D/dfa.F | 1 - source/unres/src-HCD-5D/elecont.f | 34 +- source/unres/src-HCD-5D/energy_p_new-sep_barrier.F | 124 ++-- source/unres/src-HCD-5D/energy_p_new_barrier.F | 73 ++- source/unres/src-HCD-5D/energy_split-sep.F | 76 ++- source/unres/src-HCD-5D/gen_rand_conf.F | 223 ++++++- source/unres/src-HCD-5D/gen_rand_conf_mchain.F | 424 +++++++++++++ source/unres/src-HCD-5D/geomout.F | 16 +- source/unres/src-HCD-5D/gradient_p.F | 1 + source/unres/src-HCD-5D/initialize_p.F | 182 ++++-- source/unres/src-HCD-5D/lagrangian_lesyng.F | 23 +- source/unres/src-HCD-5D/make_xx_list.F | 641 ++++++++++++++++---- source/unres/src-HCD-5D/minimize_p.F | 6 +- source/unres/src-HCD-5D/orig_frame_chain.F | 85 +++ source/unres/src-HCD-5D/readpdb-mult.F | 52 +- source/unres/src-HCD-5D/readrtns_CSA.F | 355 +++++++---- source/unres/src-HCD-5D/refsys.f | 8 +- source/unres/src-HCD-5D/stochfric.F | 25 +- source/unres/src-HCD-5D/timing.F | 37 ++ source/unres/src-HCD-5D/unres.F | 33 +- 33 files changed, 2101 insertions(+), 491 deletions(-) create mode 100644 source/unres/src-HCD-5D/gen_rand_conf_mchain.F create mode 100644 source/unres/src-HCD-5D/orig_frame_chain.F diff --git a/source/unres/src-HCD-5D/COMMON.CHAIN b/source/unres/src-HCD-5D/COMMON.CHAIN index da83764..d666fd9 100644 --- a/source/unres/src-HCD-5D/COMMON.CHAIN +++ b/source/unres/src-HCD-5D/COMMON.CHAIN @@ -1,7 +1,8 @@ integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc, & nres0,nstart_seq,nchain,chain_length,chain_border,iprzes, - & chain_border1,ireschain,tabpermchain,npermchain,afmend,afmbeg, - & nres_chomo,nmodel_start + & chain_border1,ireschain,tabpermchain,npermchain,nequiv, + & nchain_group,iequiv,mapchain,afmend,afmbeg, + & nres_chomo,nmodel_start,nran_start double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, & prod,rt,dc_work,cref,crefjlee,dc_norm2,velAFMconst, & totTafm,chomo @@ -9,14 +10,15 @@ common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), & dc_norm2(3,0:maxres2), - & dc_work(MAXRES6),nres,nres0 + & dc_work(MAXRES6),nres,nres0,nran_start common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres), & rt(3,3,maxres) common /refstruct/ cref(3,maxres2+2), & crefjlee(3,maxres2+2), & nsup,nstart_sup,nstart_seq,iprzes, & chain_length(maxchain),npermchain,ireschain(maxres), - & tabpermchain(maxchain,maxperm), + & tabpermchain(maxchain,maxperm),nchain_group, + & iequiv(maxchain,maxres),nequiv(maxchain),mapchain(maxchain), & chain_border(2,maxchain),chain_border1(2,maxchain),nchain common /from_zscore/ nz_start,nz_end,iz_sc double precision boxxsize,boxysize,boxzsize,enecut,sscut, diff --git a/source/unres/src-HCD-5D/COMMON.CONTROL b/source/unres/src-HCD-5D/COMMON.CONTROL index da8581a..32ecc93 100644 --- a/source/unres/src-HCD-5D/COMMON.CONTROL +++ b/source/unres/src-HCD-5D/COMMON.CONTROL @@ -2,7 +2,8 @@ ! and output level. !... energy_dec = .true. means print energy decomposition matrix integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad, - & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide, + & inprint,i2ndstr,mucadyn,constr_dist,symetr,npermut, + & AFMlog,selfguide, & shield_mode,tor_mode,tubelog,constr_homology,homol_nset, & iprint !... minim = .true. means DO minimization. @@ -20,6 +21,7 @@ & unres_pdb,out_cart,out_int,vdisulf,searchsc,lmuca,dccart,mucadyn, & extconf,out1file,gmatout,selfguide,AFMlog,shield_mode,tor_mode, & tubelog,constr_dist,gnorm_check,gradout,split_ene, - & with_theta_constr,with_dihed_constr,symetr,usampl,loc_qlike, + & with_theta_constr,with_dihed_constr,symetr,npermut,usampl,loc_ + & qlike, & adaptive,constr_homology,homol_nset,read2sigma,start_from_model, & read_homol_frag,out_template_coord,out_template_restr diff --git a/source/unres/src-HCD-5D/COMMON.INTERACT b/source/unres/src-HCD-5D/COMMON.INTERACT index 14416ad..9b023e5 100644 --- a/source/unres/src-HCD-5D/COMMON.INTERACT +++ b/source/unres/src-HCD-5D/COMMON.INTERACT @@ -22,22 +22,45 @@ & iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,iatel_e_vdw, & iatscp_s,iatscp_e,ispp,iscp C 3/26/20 Interaction lists +C 10/26/20 Upgraded to interaction lists for RESPA integer newcontlisti(maxint_res*maxres), + & newcontlisti_long(maxint_res*maxres), + & newcontlisti_short(maxint_res*maxres), & newcontlistj(maxint_res*maxres), + & newcontlistj_long(maxint_res*maxres), + & newcontlistj_short(maxint_res*maxres), & newcontlistppi(maxint_res*maxres), & newcontlistppj(maxint_res*maxres), - & newcontlistpp_vdwi(maxint_res*maxres), - & newcontlistpp_vdwj(maxint_res*maxres), + & newcontlistpp_vdwi_short(maxint_res*maxres), + & newcontlistpp_vdwj_short(maxint_res*maxres), & newcontlistscpi(2*maxint_res*maxres), + & newcontlistscpi_long(2*maxint_res*maxres), + & newcontlistscpi_short(2*maxint_res*maxres), & newcontlistscpj(2*maxint_res*maxres), + & newcontlistscpj_long(2*maxint_res*maxres), + & newcontlistscpj_short(2*maxint_res*maxres), & g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_end, - & g_listpp_vdw_start,g_listpp_vdw_end,g_listscp_start,g_listscp_end - common /interact_list/newcontlisti,newcontlistj,g_listscsc_start, - & g_listscsc_end,newcontlistppi,newcontlistppj,g_listpp_start, - & g_listpp_end,newcontlistpp_vdwi,newcontlistpp_vdwj, - & g_listpp_vdw_start,g_listpp_vdw_end, - & newcontlistscpi,newcontlistscpj,g_listscp_start, - & g_listscp_end + & g_listscp_start,g_listscp_end, + & g_listscsc_start_long,g_listscsc_end_long, + & g_listscp_start_long,g_listscp_end_long, + & g_listscsc_start_short,g_listscsc_end_short, + & g_listpp_vdw_start_short,g_listpp_vdw_end_short, + & g_listscp_start_short,g_listscp_end_short + common /interact_list/ + & newcontlisti,newcontlisti_long,newcontlisti_short, + & newcontlistj,newcontlistj_long,newcontlistj_short, + & g_listscsc_start,g_listscsc_start_long,g_listscsc_start_short, + & g_listscsc_end,g_listscsc_end_long,g_listscsc_end_short, + & newcontlistppi,newcontlistppj,g_listpp_start, + & g_listpp_end, + & newcontlistpp_vdwi_short, + & newcontlistpp_vdwj_short, + & g_listpp_vdw_start_short, + & g_listpp_vdw_end_short, + & newcontlistscpi,newcontlistscpi_long,newcontlistscpi_short, + & newcontlistscpj,newcontlistscpj_long,newcontlistscpj_short, + & g_listscp_start,g_listscp_start_long,g_listscp_start_short, + & g_listscp_end,g_listscp_end_long,g_listscp_end_short C 12/1/95 Array EPS included in the COMMON block. double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp, & sigma0,sigii, diff --git a/source/unres/src-HCD-5D/COMMON.MD b/source/unres/src-HCD-5D/COMMON.MD index cea18eb..c69e8e4 100644 --- a/source/unres/src-HCD-5D/COMMON.MD +++ b/source/unres/src-HCD-5D/COMMON.MD @@ -3,13 +3,13 @@ & dvmax,damax,edriftmax integer n_timestep,ntwx,ntwe,lang,count_reset_moment, & count_reset_vel,ntime_split,ntime_split0, - & maxtime_split,itime_mat,imatupdate + & maxtime_split,itime_mat,imatupdate,irest_freq logical large,print_compon,tbf,rest,reset_moment,reset_vel, & rattle,mdpdb,RESPA,preminim common /mdpar/ v_ini,d_time,d_time0,t_bath, & tau_bath,dvmax,damax,n_timestep,mdpdb, & ntime_split,ntime_split0,maxtime_split,itime_mat,imatupdate, - & ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim, + & ntwx,ntwe,irest_freq,lang,large,print_compon,tbf,rest,preminim, & reset_moment,reset_vel,count_reset_moment,count_reset_vel, & rattle,RESPA ! Basic quantities diff --git a/source/unres/src-HCD-5D/COMMON.TIME1 b/source/unres/src-HCD-5D/COMMON.TIME1 index fc1e7d5..24d82ee 100644 --- a/source/unres/src-HCD-5D/COMMON.TIME1 +++ b/source/unres/src-HCD-5D/COMMON.TIME1 @@ -15,7 +15,10 @@ c FOUND_NAN - set by calcf to stop sumsl via stopx & time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric, & time_scatter_fmat,time_scatter_ginv, & time_fsample,time_scatter_fmatmult,time_scatter_ginvmult, - & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS + & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list, + & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short, + & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long, + & time_escpcalc,time_escpsetup common /timing/ t_init,t_MDsetup,t_langsetup, & t_MD,t_enegrad,t_sdsetup,time_bcast,time_reduce,time_gather, & time_sendrecv,time_scatter,time_barrier_e,time_barrier_g, @@ -25,4 +28,7 @@ c FOUND_NAN - set by calcf to stop sumsl via stopx & time_vec,time_mat,time_ginvmult,time_fricmatmult,time_fric, & time_fsample,time_scatter_fmatmult,time_scatter_ginvmult, & time_scatter_fmat,time_scatter_ginv, - & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS + & time_stoch,t_eshort,t_elong,t_etotal,time_SAXS,time_list, + & time_evdw,time_eelec,time_escp,time_evdw_short,time_eelec_short, + & time_escp_short,time_evdw_long,time_eelec_long,time_escp_long, + & time_escpcalc,time_escpsetup diff --git a/source/unres/src-HCD-5D/COMMON.TORCNSTR b/source/unres/src-HCD-5D/COMMON.TORCNSTR index 9476b50..5219a92 100644 --- a/source/unres/src-HCD-5D/COMMON.TORCNSTR +++ b/source/unres/src-HCD-5D/COMMON.TORCNSTR @@ -1,4 +1,5 @@ - integer ndih_constr,idih_constr(maxdih_constr),ntheta_constr, + integer ndih_constr,idih_constr(maxdih_constr), + & iconstr_dih(maxres),ntheta_constr, & itheta_constr(maxdih_constr) integer ndih_nconstr,idih_nconstr(maxdih_constr) integer idihconstr_start,idihconstr_end,ithetaconstr_start, @@ -10,7 +11,7 @@ & vpsipred(3,maxdih_constr),sdihed(2,maxdih_constr),wdihc common /torcnstr/ phi0,drange,ftors,theta_constr0,theta_drange, & for_thet_constr,vpsipred,sdihed,wdihc, - & ndih_constr,idih_constr, + & ndih_constr,idih_constr,iconstr_dih, & ndih_nconstr,idih_nconstr,idihconstr_start,idihconstr_end, & ntheta_constr,itheta_constr,ithetaconstr_start, & ithetaconstr_end,raw_psipred diff --git a/source/unres/src-HCD-5D/DIMENSIONS b/source/unres/src-HCD-5D/DIMENSIONS index 4236776..3c2c924 100644 --- a/source/unres/src-HCD-5D/DIMENSIONS +++ b/source/unres/src-HCD-5D/DIMENSIONS @@ -16,7 +16,8 @@ C Max. number of coarse-grain processors parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=20000) +c parameter (maxres=70000) + parameter (maxres=1000) C Max. number of AA residues per chain integer maxres_chain parameter (maxres_chain=1200) @@ -30,9 +31,10 @@ C Appr. max. number of interaction sites & mmaxres2_chain=maxres2_chain*(maxres2_chain+1)/2) C Max number of symetric chains integer maxchain - parameter (maxchain=50) + parameter (maxchain=250) integer maxperm - parameter (maxperm=5040) +c parameter (maxperm=5040) + parameter (maxperm=3628800) C Max. number of variables integer maxvar parameter (maxvar=6*maxres) @@ -42,8 +44,8 @@ C Max. number of groups of interactions that a given SC is involved in C Max. number of derivatives of virtual-bond and side-chain vectors in theta C or phi. integer maxdim -c parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2) - parameter (maxdim=(maxres-1)*(maxres-2)/2) + parameter (maxdim=(maxres_chain-1)*(maxres_chain-2)/2) +c parameter (maxdim=(maxres-1)*(maxres-2)/2) C Max. number of SC contacts integer maxcont parameter (maxcont=12*maxres) @@ -60,8 +62,7 @@ C include in template-based/contact distance restraints. parameter (maxcont_res=200) C Max. number of distance/contact-distance restraints integer maxdim_cont -c parameter (maxdim_cont=maxres*maxcont_res) - parameter (maxdim_cont=maxres*1000) + parameter (maxdim_cont=maxres*maxcont_res) C Number of AA types (at present only natural AA's will be handled integer ntyp,ntyp1 parameter (ntyp=24,ntyp1=ntyp+1) diff --git a/source/unres/src-HCD-5D/MD_A-MTS.F b/source/unres/src-HCD-5D/MD_A-MTS.F index e504cbd..8cebc8f 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -65,6 +65,7 @@ c t_enegrad=0.0d0 t_sdsetup=0.0d0 write (iout,'(20(1h=),a20,20(1h=))') "MD calculation started" + write (iout,*) "Restart frequency:",irest_freq #ifdef MPI tt0=MPI_Wtime() #else @@ -261,7 +262,7 @@ C call check_ecartint call cartout(totT) endif endif - if (rstcount.eq.1000.or.itime.eq.n_timestep) then + if (rstcount.eq.irest_freq.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath do i=0,2*nres-1 @@ -290,7 +291,7 @@ C call check_ecartint & ' End of MD calculation ' #ifdef TIMING_ENE write (iout,*) "time for etotal",t_etotal," elong",t_elong, - & " eshort",t_eshort + & " eshort",t_eshort," list",time_list write (iout,*) "time_fric",time_fric," time_stoch",time_stoch, & " time_fricmatmult",time_fricmatmult," time_fsample ", & time_fsample @@ -429,7 +430,7 @@ c Calculate energy and forces call zerograd call etotal(potEcomp) ! AL 4/17/17: Reduce the steps if NaNs occurred. - if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then + if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0))) then d_time=d_time/2 cycle endif @@ -1672,6 +1673,7 @@ c Set up the initial conditions of a MD simulation double precision etot logical fail integer i_start_models(0:nodes-1) + double precision potEcomp_long(0:Max_Ene) write (iout,*) "init_MD INDPDB",indpdb d_time0=d_time c write(iout,*) "d_time", d_time @@ -1810,7 +1812,7 @@ c rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then write (iout,*) "Initial velocities" do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + write (iout,'(i7,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo endif @@ -2080,6 +2082,7 @@ C 8/22/17 AL Minimize initial structure call zerograd call etotal(potEcomp) if (large) call enerprint(potEcomp) + if (RESPA) call etotal_long(potEcomp_long) #ifdef TIMING_ENE #ifdef MPI t_etotal=t_etotal+MPI_Wtime()-tt0 @@ -2092,7 +2095,7 @@ c write (iout,*) "PotE-homology",potE-potEcomp(27) call cartgrad call lagrangian call max_accel - if (amax*d_time .gt. dvmax) then + if (amax*d_time .gt. dvmax .and. .not. respa) then d_time=d_time*dvmax/amax if(me.eq.king.or..not.out1file) write (iout,*) & "Time step reduced to",d_time, diff --git a/source/unres/src-HCD-5D/cartprint.f b/source/unres/src-HCD-5D/cartprint.f index 339a89d..8a747b4 100644 --- a/source/unres/src-HCD-5D/cartprint.f +++ b/source/unres/src-HCD-5D/cartprint.f @@ -15,6 +15,6 @@ & ' centroid coordinates'/ 1 ' ', 7X,'X',11X,'Y',11X,'Z', & 10X,'X',11X,'Y',11X,'Z') - 110 format (a,'(',i4,')',6f12.5) + 110 format (a,'(',i5,')',6f12.5) return end diff --git a/source/unres/src-HCD-5D/chain_symmetry.F b/source/unres/src-HCD-5D/chain_symmetry.F index 8c36855..d10c33f 100644 --- a/source/unres/src-HCD-5D/chain_symmetry.F +++ b/source/unres/src-HCD-5D/chain_symmetry.F @@ -1,5 +1,6 @@ subroutine chain_symmetry(nchain,nres,itype,chain_border, - & chain_length,npermchain,tabpermchain) + & chain_length,npermchain,tabpermchain,nchain_group,nequiv, + & iequiv,chaingroup) c c Determine chain symmetry. nperm is the number of permutations and c tabperchain contains the allowed permutations of the chains. @@ -12,13 +13,17 @@ c & chain_length(nchain),itemp(maxchain), & npermchain,tabpermchain(maxchain,maxperm), & tabperm(maxchain,maxperm),mapchain(maxchain), - & iequiv(maxchain,maxres),iflag(maxres) + & chaingroup(maxchain),iequiv(maxchain,maxres),iflag(maxres) integer i,j,k,l,ii,nchain_group,nequiv(maxchain),iieq, & nperm,npermc,ind + logical lprn /.false./ if (nchain.eq.1) then npermchain=1 tabpermchain(1,1)=1 c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1) + nchain_group=1 + iequiv(1,1)=1 + chaingroup(1)=1 return endif c @@ -73,12 +78,24 @@ c k=k+1 do j=1,nequiv(i) ind=ind+1 mapchain(ind)=iequiv(j,i) + chaingroup(ind)=i enddo enddo write (iout,*) "mapchain" do i=1,nchain write (iout,*) i,mapchain(i) enddo + write (iout,*) "chaingroup" + do i=1,nchain + write (iout,*) i,chaingroup(i) + enddo + if (npermut.eq.0) then + npermchain=1 + do i=1,nchain + tabpermchain(i,1)=i + enddo + return + endif ii=0 do i=1,nchain_group call permut(nequiv(i),nperm,tabperm) @@ -117,10 +134,12 @@ c k=k+1 enddo enddo write(iout,*) "Number of chain permutations",npermchain + if (lprn) then write(iout,*) "Permutations" do i=1,npermchain write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain) enddo + endif return end c--------------------------------------------------------------------- diff --git a/source/unres/src-HCD-5D/chainbuild.F b/source/unres/src-HCD-5D/chainbuild.F index 7902f15..a60e2bd 100644 --- a/source/unres/src-HCD-5D/chainbuild.F +++ b/source/unres/src-HCD-5D/chainbuild.F @@ -108,11 +108,11 @@ C include 'COMMON.VAR' cost=dcos(theta(3)) sint=dsin(theta(3)) - t(1,1,1)=-cost + t(1,1,1)=cost t(1,2,1)=-sint t(1,3,1)= 0.0D0 - t(2,1,1)=-sint - t(2,2,1)= cost + t(2,1,1)=sint + t(2,2,1)=cost t(2,3,1)= 0.0D0 t(3,1,1)= 0.0D0 t(3,2,1)= 0.0D0 @@ -170,7 +170,7 @@ C C Locate CA(i) and SC(i-1) C implicit none - integer i,j + integer i,j,k double precision theti,phii,cost,sint,cosphi,sinphi include 'DIMENSIONS' include 'COMMON.CHAIN' @@ -196,14 +196,15 @@ C sint=dsin(theti) cosphi=dcos(phii) sinphi=dsin(phii) +c write (iout,*) "locate_next_res i",i * Define the matrices of the rotation about the virtual-bond valence angles * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this * program), R(i,j,k), and, the cumulative matrices of rotation RT t(1,1,i-2)=-cost - t(1,2,i-2)=-sint + t(1,2,i-2)= sint t(1,3,i-2)= 0.0D0 t(2,1,i-2)=-sint - t(2,2,i-2)= cost + t(2,2,i-2)=-cost t(2,3,i-2)= 0.0D0 t(3,1,i-2)= 0.0D0 t(3,2,i-2)= 0.0D0 @@ -212,26 +213,36 @@ C r(1,2,i-2)= 0.0D0 r(1,3,i-2)= 0.0D0 r(2,1,i-2)= 0.0D0 - r(2,2,i-2)=-cosphi - r(2,3,i-2)= sinphi + r(2,2,i-2)= cosphi + r(2,3,i-2)=-sinphi r(3,1,i-2)= 0.0D0 r(3,2,i-2)= sinphi r(3,3,i-2)= cosphi rt(1,1,i-2)=-cost - rt(1,2,i-2)=-sint + rt(1,2,i-2)= sint rt(1,3,i-2)=0.0D0 - rt(2,1,i-2)=sint*cosphi + rt(2,1,i-2)=-sint*cosphi rt(2,2,i-2)=-cost*cosphi - rt(2,3,i-2)=sinphi + rt(2,3,i-2)=-sinphi rt(3,1,i-2)=-sint*sinphi - rt(3,2,i-2)=cost*sinphi - rt(3,3,i-2)=cosphi + rt(3,2,i-2)=-cost*sinphi + rt(3,3,i-2)= cosphi call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1)) +c write (iout,*) "prod",i-2 +c do j=1,3 +c write (iout,*) (prod(j,k,i-2),k=1,3) +c enddo +c write (iout,*) "prod",i-1 +c do j=1,3 +c write (iout,*) (prod(j,k,i-1),k=1,3) +c enddo do j=1,3 dc_norm(j,i-1)=prod(j,1,i-1) dc(j,i-1)=vbld(i)*prod(j,1,i-1) c(j,i)=c(j,i-1)+dc(j,i-1) enddo +c write (iout,*) "dc",i-1,(dc(j,i-1),j=1,3) +c write (iout,*) "c",i,(dc(j,i),j=1,3) cd print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3) C C Now calculate the coordinates of SC(i-1) diff --git a/source/unres/src-HCD-5D/checkder_p.F b/source/unres/src-HCD-5D/checkder_p.F index e1448db..d0032da 100644 --- a/source/unres/src-HCD-5D/checkder_p.F +++ b/source/unres/src-HCD-5D/checkder_p.F @@ -364,6 +364,7 @@ c------------------------------------------------------------------------- double precision dnorm1,dnorm2,be double precision time00 double precision dist,alpha,beta + double precision time01 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates' #ifdef TIMING time01=MPI_Wtime() diff --git a/source/unres/src-HCD-5D/contact.f b/source/unres/src-HCD-5D/contact.f index f986380..ead0332 100644 --- a/source/unres/src-HCD-5D/contact.f +++ b/source/unres/src-HCD-5D/contact.f @@ -137,7 +137,7 @@ c print *,'nnt=',nnt,' nct=',nct i2=icont(2,i) it1=itype(i1) it2=itype(i2) - write (iout,'(i5,2x,a,i5,2x,a,i4)') + write (iout,'(i10,2x,a,i7,2x,a,i7)') & i,restyp(it1),i1,restyp(it2),i2 enddo endif @@ -184,8 +184,8 @@ c enddo ii1=iharp(3,i) jj1=iharp(4,i) write (iout,*) - write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=i1,ii1) - write (iout,'(20(a,i5,1x))') (restyp(itype(k)),k,k=j1,jj1,-1) + write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=i1,ii1) + write (iout,'(20(a,i7,1x))') (restyp(itype(k)),k,k=j1,jj1,-1) c do k=jj1,j1,-1 c write (iout,'(a,i3,$)') restyp(itype(k)),k c enddo diff --git a/source/unres/src-HCD-5D/dfa.F b/source/unres/src-HCD-5D/dfa.F index 1af0b44..3982b6d 100644 --- a/source/unres/src-HCD-5D/dfa.F +++ b/source/unres/src-HCD-5D/dfa.F @@ -368,7 +368,6 @@ C END OF BETA RESTRAINT edfadis=0 gdfad=0.0d0 -c write (2,*) "edfad",idfadis_start,idfadis_end do i=idfadis_start,idfadis_end iatm1=idislis(1,i)+ishiftca diff --git a/source/unres/src-HCD-5D/elecont.f b/source/unres/src-HCD-5D/elecont.f index bf9056a..7c024ea 100644 --- a/source/unres/src-HCD-5D/elecont.f +++ b/source/unres/src-HCD-5D/elecont.f @@ -9,7 +9,7 @@ include 'COMMON.NAMES' logical lprint double precision elpp_6(2,2),elpp_3(2,2),ael6_(2,2),ael3_(2,2) - double precision app_(2,2),bpp_(2,2),rpp_(2,2) + double precision app_(2,2),bpp_(2,2),epp_(2,2),rpp_(2,2) integer ncont,icont(2,maxcont) double precision econt(maxcont) double precision boxshift @@ -20,7 +20,7 @@ * * as of 7/06/91. * -c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ +c data epp_ / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ data rpp_ / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/ data elpp_6 /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/ data elpp_3 / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/ @@ -55,7 +55,6 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ zmedi=zi+0.5*dzi call to_box(xmedi,ymedi,zmedi) c write (iout,*) "i",xmedi,ymedi,zmedi -c write (iout,*) "i",xmedi,ymedi,zmedi do 4 j=i+2,nct-1 c write (iout,*) "i",i," j",j if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) goto 4 @@ -74,11 +73,12 @@ c write (iout,*) "i",i," j",j xj=c(1,j)+0.5*dxj yj=c(2,j)+0.5*dyj zj=c(3,j)+0.5*dzj -c write (iout,*) "j",xj,yj,zj call to_box(xj,yj,zj) - xj=boxshift(xj-xi,boxxsize) - yj=boxshift(yj-yi,boxysize) - zj=boxshift(zj-zi,boxzsize) +c write (iout,*) "j",xj,yj,zj + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) +c write (iout,*) "j",xj,yj,zj rij=xj*xj+yj*yj+zj*zj rrmij=1.0/(xj*xj+yj*yj+zj*zj) rmij=sqrt(rrmij) @@ -105,7 +105,7 @@ c write (iout,*) "j",xj,yj,zj econt(ncont)=eesij endif ees=ees+eesij -c write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij +c write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij 4 continue 1 continue if (lprint) then @@ -118,7 +118,7 @@ c write (iout,*) "i"," j",j," rij",dsqrt(rij)," eesij",eesij i2=icont(2,i) it1=itype(i1) it2=itype(i2) - write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') + write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif @@ -204,7 +204,7 @@ c enddo i2=icont(2,i) it1=itype(i1) it2=itype(i2) - write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)') + write (iout,'(i7,2x,a,i6,2x,a,i6,f10.5)') & i,restyp(it1),i1,restyp(it2),i2,econt(i) enddo endif @@ -267,7 +267,7 @@ cd write (iout,*) i1,j1,not_done ii1=max0(ii1-1,1) jj1=max0(jj1-1,1) nbeta=nbeta+1 - if(lprint)write(iout,'(a,i3,4i4)')'parallel beta', + if(lprint)write(iout,'(a,i7,4i6)')'parallel beta', & nbeta,ii1,i1,jj1,j1 nbfrag=nbfrag+1 @@ -442,29 +442,29 @@ cd write (iout,*) i1,j1,not_done if (lprint) then - write (iout,'(a,i3,4i4)')'antiparallel beta', + write (iout,'(a,i3,4i6)')'antiparallel beta', & nbeta,ii1-1,i1,jj1,j1-1 nstrand=nstrand+1 if (nstrand.le.9) then - write(12,'(a18,i1,a9,i3,a2,i3,a1)') + write(12,'(a18,i1,a9,i6,a2,i6,a1)') & "DefPropRes 'strand",nstrand, & "' 'num = ",ii1-2,"..",i1-1,"'" else - write(12,'(a18,i2,a9,i3,a2,i3,a1)') + write(12,'(a18,i2,a9,i6,a2,i6,a1)') & "DefPropRes 'strand",nstrand, & "' 'num = ",ii1-2,"..",i1-1,"'" endif nstrand=nstrand+1 if (nstrand.le.9) then - write(12,'(a18,i1,a9,i3,a2,i3,a1)') + write(12,'(a18,i1,a9,i6,a2,i6,a1)') & "DefPropRes 'strand",nstrand, & "' 'num = ",j1-2,"..",jj1-1,"'" else - write(12,'(a18,i2,a9,i3,a2,i3,a1)') + write(12,'(a18,i2,a9,i6,a2,i6,a1)') & "DefPropRes 'strand",nstrand, & "' 'num = ",j1-2,"..",jj1-1,"'" endif - write(12,'(a8,4i4)') + write(12,'(a8,4i6)') & "SetNeigh",ii1-2,i1-1,jj1-1,j1-2 endif endif diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F index c4e54bc..f92aebb 100644 --- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F @@ -94,9 +94,9 @@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon gg_lipi=0.0d0 gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -219,9 +219,9 @@ c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon gg_lipi=0.0d0 gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -340,9 +340,9 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon gg_lipi=0.0d0 gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -462,9 +462,9 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon gg_lipi=0.0d0 gg_lipj=0.0d0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -586,9 +586,9 @@ c else c endif ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -734,9 +734,9 @@ c else c endif ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -872,7 +872,7 @@ C & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision subchap,sss1,sssgrad1 - double precision boxshift + double precision boxshift,rij1 evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -882,9 +882,12 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + if (energy_dec) + & write(2,*) "g_listscsc_start_long,g_listscsc_end_long", + & g_listscsc_start_long,g_listscsc_end_long + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -940,9 +943,13 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) - sss1=sscale(1.0d0/rij,r_cut_int) + rij1=1.0d0/rij +c sss1=sscale(1.0d0/rij,r_cut_int) + sss1=sscale(rij1,r_cut_int) if (sss1.eq.0.0d0) cycle - sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) + rij1=rij1/sigmaii(itypi,itypj) + sss=sscale(rij1,r_cut_respa) +c sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) if (sss.lt.1.0d0) then C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -1026,6 +1033,7 @@ C This subroutine calculates the interaction energy of nonbonded side chains C assuming the Gay-Berne potential of interaction. C implicit none + include 'mpif.h' include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -1038,6 +1046,7 @@ C include 'COMMON.CALC' include 'COMMON.CONTROL' include "COMMON.SPLITELE" + include 'COMMON.TIME1' logical lprn double precision evdw integer itypi,itypj,itypi1,iint,ind,ikont @@ -1046,6 +1055,8 @@ C & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip double precision dist,sscale,sscagrad,sscagradlip,sscalelip double precision boxshift + double precision time01 +c time01=MPI_Wtime() evdw=0.0D0 ccccc energy_dec=.false. c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -1055,9 +1066,12 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.false. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + if (energy_dec) + & write(2,*) "g_listscsc_start_short,g_listscsc_end_short", + & g_listscsc_start_short,g_listscsc_end_short + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1119,8 +1133,8 @@ c & (2.0d0-sslipi-sslipj)/2.0d0 rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa) - sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) if (sss.gt.0.0d0) then + sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa) C Calculate angle-dependent terms of energy and contributions to their C derivatives. @@ -1189,6 +1203,7 @@ C Calculate angular part of the gradient. c enddo ! j c enddo ! iint enddo ! i +c time_evdw_short=time_evdw_short+MPI_Wtime()-time01 c write (iout,*) "Number of loop steps in EGB:",ind cccc energy_dec=.false. return @@ -1230,9 +1245,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_long,g_listscsc_end_long + i=newcontlisti_long(ikont) + j=newcontlistj_long(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1392,9 +1407,9 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon c if (icall.eq.0) lprn=.true. ind=0 c do i=iatsc_s,iatsc_e - do ikont=g_listscsc_start,g_listscsc_end - i=newcontlisti(ikont) - j=newcontlistj(ikont) + do ikont=g_listscsc_start_short,g_listscsc_end_short + i=newcontlisti_short(ikont) + j=newcontlistj_short(ikont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1742,6 +1757,9 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c c do i=iatel_s,iatel_e + if (energy_dec) + & write(iout,*) "g_listpp_start,g_listpp_end", + & g_listpp_start,g_listpp_end do ikont=g_listpp_start,g_listpp_end i=newcontlistppi(ikont) j=newcontlistppj(ikont) @@ -2812,9 +2830,12 @@ c write (iout,*) "iatel_s_vdw",iatel_s_vdw, c & " iatel_e_vdw",iatel_e_vdw c call flush(iout) c do i=iatel_s_vdw,iatel_e_vdw - do ikont=g_listpp_vdw_start,g_listpp_vdw_end - i=newcontlistpp_vdwi(ikont) - j=newcontlistpp_vdwj(ikont) + if (energy_dec) + & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short", + & g_listpp_vdw_start_short,g_listpp_vdw_end_short + do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short + i=newcontlistpp_vdwi_short(ikont) + j=newcontlistpp_vdwj_short(ikont) if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) @@ -2939,9 +2960,12 @@ c if (lprint_short) c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e c do i=iatscp_s,iatscp_e - do ikont=g_listscp_start,g_listscp_end - i=newcontlistscpi(ikont) - j=newcontlistscpj(ikont) + if (energy_dec) + & write(iout,*)"g_listscp_start_long,g_listscp_end_long", + & g_listscp_start_long,g_listscp_end_long + do ikont=g_listscp_start_long,g_listscp_end_long + i=newcontlistscpi_long(ikont) + j=newcontlistscpj_long(ikont) if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) @@ -3065,14 +3089,20 @@ C double precision ggg(3) double precision sscale,sscagrad double precision boxshift + integer ikont evdw2=0.0D0 evdw2_14=0.0d0 cd print '(a)','Enter ESCP' c if (lprint_short) c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s, c & ' iatscp_e=',iatscp_e - if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb - do i=iatscp_s,iatscp_e +c if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb + if (energy_dec) + & write(iout,*) "g_listscp_start_short,g_listscp_end_short", + & g_listscp_start_short,g_listscp_end_short + do ikont=g_listscp_start_short,g_listscp_end_short + i=newcontlistscpi_short(ikont) + j=newcontlistscpj_short(ikont) if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) @@ -3083,9 +3113,9 @@ c & ' iatscp_e=',iatscp_e c if (lprint_short) c & write (iout,*) "i",i," itype",itype(i),itype(i+1), c & " nscp_gr",nscp_gr(i) - do iint=1,nscp_gr(i) - - do j=iscpstart(i,iint),iscpend(i,iint) +c do iint=1,nscp_gr(i) +c +c do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) c if (lprint_short) c & write (iout,*) "j",j," itypj",itypj @@ -3149,9 +3179,9 @@ c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo endif - enddo +c enddo - enddo ! iint +c enddo ! iint enddo ! i do i=1,nct do j=1,3 diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index 6d6a817..2d94dc0 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -37,6 +37,7 @@ c include 'COMMON.MD' & eliptran,Eafmforce,Etube, & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet integer n_corr,n_corr1 + double precision time01 #ifdef MPI c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, c & " nfgtasks",nfgtasks @@ -117,6 +118,9 @@ c call chainbuild_cart endif c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate if (mod(itime_mat,imatupdate).eq.0) then +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif call make_SCp_inter_list c write (iout,*) "Finished make_SCp_inter_list" c call flush(iout) @@ -126,9 +130,12 @@ c call flush(iout) call make_pp_inter_list c write (iout,*) "Finished make_pp_inter_list" c call flush(iout) - call make_pp_vdw_inter_list +c call make_pp_vdw_inter_list c write (iout,*) "Finished make_pp_vdw_inter_list" c call flush(iout) +#ifdef TIMING_ENE + time_list=time_list+MPI_Wtime()-time01 +#endif endif c print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct @@ -151,6 +158,9 @@ C C Compute the side-chain and electrostatic interaction energy C C print *,ipot +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj(evdw) @@ -175,6 +185,9 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef TIMING_ENE + time_evdw=time_evdw+MPI_Wtime()-time01 +#endif #ifdef DFA C BARTEK for dfa test! c print *,"Processors",MyRank," wdfa",wdfa_dist @@ -216,6 +229,9 @@ c print *,"Processor",myrank," computed USCSC" #ifdef TIMING time_vec=time_vec+MPI_Wtime()-time01 #endif +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif C Introduction of shielding effect first for each peptide group C the shielding factor is set this factor is describing how each C peptide group is shielded by side-chains @@ -252,6 +268,9 @@ c print *,"Processor",myrank," left VEC_AND_DERIV" c call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, c & eello_turn4) endif +#ifdef TIMING_ENE + time_eelec=time_eelec+MPI_Wtime()-time01 +#endif c#ifdef TIMING c time_enecalc=time_enecalc+MPI_Wtime()-time00 c#endif @@ -260,6 +279,9 @@ C C Calculate excluded-volume interaction energy between peptide groups C and side chains. C +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif if (ipot.lt.6) then if(wscp.gt.0d0) then call escp(evdw2,evdw2_14) @@ -271,6 +293,9 @@ C c write (iout,*) "Soft-sphere SCP potential" call escp_soft_sphere(evdw2,evdw2_14) endif +#ifdef TIMING_ENE + time_escp=time_escp+MPI_Wtime()-time01 +#endif c c Calculate the bond-stretching energy c @@ -839,12 +864,12 @@ c do i=nnt,nres gradbufc(k,i)=0.0d0 enddo enddo -#ifdef DEBUG - write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end - write (iout,*) (i," jgrad_start",jgrad_start(i), - & " jgrad_end ",jgrad_end(i), - & i=igrad_start,igrad_end) -#endif +c#ifdef DEBUG +c write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end +c write (iout,*) (i," jgrad_start",jgrad_start(i), +c & " jgrad_end ",jgrad_end(i), +c & i=igrad_start,igrad_end) +c#endif c c Obsolete and inefficient code; we can make the effort O(n) and, therefore, c do not parallelize this part. @@ -1510,6 +1535,7 @@ C double precision sscale,sscagrad,sscagradlip,sscalelip double precision gg_lipi(3),gg_lipj(3) double precision boxshift + external boxshift c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 gg_lipi=0.0d0 @@ -2039,6 +2065,7 @@ c do j=istart(i,iint),iend(i,iint) c write(iout,*) "PRZED ZWYKLE", evdwij call dyn_ssbond_ene(i,j,evdwij) c write(iout,*) "PO ZWYKLE", evdwij +c call flush(iout) evdw=evdw+evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') @@ -2151,7 +2178,7 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') cd & restyp(itypi),i,restyp(itypj),j, cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) -c return + return endif sigder=-sig*sigsq c--------------------------------------------------------------- @@ -3049,7 +3076,7 @@ c write(iout,*) 'b2=',(b2(k,i-2),k=1,2) #endif enddo - mu=0.0d0 + mu(:,:nres)=0.0d0 #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else @@ -5504,6 +5531,9 @@ C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C implicit none +#ifdef MPI + include 'mpif.h' +#endif include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -5515,6 +5545,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.SPLITELE' + include 'COMMON.TIME1' double precision ggg(3) integer i,iint,j,k,iteli,itypj,subchap,ikont double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, @@ -5522,6 +5553,10 @@ C double precision evdw2,evdw2_14,evdwij double precision sscale,sscagrad double precision boxshift + external boxshift,to_box +c#ifdef TIMING_ENE +c double precision time01 +c#endif evdw2=0.0D0 evdw2_14=0.0d0 c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' @@ -5533,6 +5568,9 @@ C do zshift=-1,1 if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb c do i=iatscp_s,iatscp_e do ikont=g_listscp_start,g_listscp_end +c#ifdef TIMING_ENE +c time01=MPI_Wtime() +c#endif i=newcontlistscpi(ikont) j=newcontlistscpj(ikont) if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle @@ -5540,6 +5578,7 @@ c do i=iatscp_s,iatscp_e xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) zi=0.5D0*(c(3,i)+c(3,i+1)) +!DIR$ INLINE call to_box(xi,yi,zi) c do iint=1,nscp_gr(i) @@ -5554,11 +5593,21 @@ C Uncomment following three lines for Ca-p interactions xj=c(1,j) yj=c(2,j) zj=c(3,j) +!DIR$ INLINE call to_box(xj,yj,zj) +c#ifdef TIMING_ENE +c time_escpsetup=time_escpsetup+MPI_Wtime()-time01 +c time01=MPI_Wtime() +c#endif +!DIR$ INLINE xj=boxshift(xj-xi,boxxsize) yj=boxshift(yj-yi,boxysize) zj=boxshift(zj-zi,boxzsize) c print *,xj,yj,zj,'polozenie j' +c#ifdef TIMING_ENE +c time_escpsetup=time_escpsetup+MPI_Wtime()-time01 +c time01=MPI_Wtime() +c#endif rrij=1.0D0/(xj*xj+yj*yj+zj*zj) c print *,rrij sss=sscale(1.0d0/(dsqrt(rrij)),r_cut_int) @@ -5619,6 +5668,9 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo +c#ifdef TIMING_ENE +c time_escpcalc=time_escpcalc+MPI_Wtime()-time01 +c#endif c endif !endif for sscale cutoff c enddo ! j @@ -5804,7 +5856,8 @@ C 15/02/13 CC dynamic SSbond - additional check if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) - ehpb=ehpb+2*eij +c ehpb=ehpb+2*eij + ehpb=ehpb+eij endif cd write (iout,*) "eij",eij cd & ' waga=',waga,' fac=',fac diff --git a/source/unres/src-HCD-5D/energy_split-sep.F b/source/unres/src-HCD-5D/energy_split-sep.F index 34a1bd1..51d5b2d 100644 --- a/source/unres/src-HCD-5D/energy_split-sep.F +++ b/source/unres/src-HCD-5D/energy_split-sep.F @@ -13,7 +13,7 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include "mpif.h" double precision weights_(n_ene) - double precision time00,time_Bcast,time_BcastW + double precision time00 integer ierror,ierr #endif include 'COMMON.SETUP' @@ -29,6 +29,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.QRESTR' include 'COMMON.MD' include 'COMMON.CONTROL' + include 'COMMON.TIME1' double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6, @@ -36,6 +37,9 @@ cMS$ATTRIBUTES C :: proc_proc & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet integer i,n_corr,n_corr1 c write(iout,'(a,i2)')'Calling etotal_long ipot=',ipot +#ifdef TIMING_ENE + double precision time01 +#endif if (modecalc.eq.12.or.modecalc.eq.14) then #ifdef MPI c if (fg_rank.eq.0) call int_from_cart1(.false.) @@ -138,16 +142,25 @@ c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) endif if (mod(itime_mat,imatupdate).eq.0) then - call make_SCp_inter_list - call make_SCSC_inter_list +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif + call make_SCp_inter_list_RESPA + call make_SCSC_inter_list_RESPA call make_pp_inter_list - call make_pp_vdw_inter_list + call make_pp_vdw_inter_list_RESPA +#ifdef TIMING_ENE + time_list=time_list+MPI_Wtime()-time01 +#endif endif #endif cd print *,'nnt=',nnt,' nct=',nct C C Compute the side-chain and electrostatic interaction energy C +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj_long(evdw) @@ -171,8 +184,20 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef TIMING_ENE + time_evdw_long=time_evdw_long+MPI_Wtime()-time01 +#endif +#ifdef TIMING + time01=MPI_Wtime() +#endif call vec_and_deriv +#ifdef TIMING + time_vec=time_vec+MPI_Wtime()-time01 +#endif c write (iout,*) "etotal_long: shield_mode",shield_mode +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif if (shield_mode.eq.1) then call set_shield_fac else if (shield_mode.eq.2) then @@ -204,10 +229,16 @@ c write (iout,*) "Soft-spheer ELEC potential" call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, & eello_turn4) endif +#ifdef TIMING_ENE + time_eelec_long=time_eelec_long+MPI_Wtime()-time01 +#endif C C Calculate excluded-volume interaction energy between peptide groups C and side chains. C +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif if (ipot.lt.6) then if(wscp.gt.0d0) then call escp_long(evdw2,evdw2_14) @@ -218,6 +249,9 @@ C else call escp_soft_sphere(evdw2,evdw2_14) endif +#ifdef TIMING_ENE + time_escp_long=time_escp_long+MPI_Wtime()-time01 +#endif #ifdef FOURBODY C C 12/1/95 Multi-body terms @@ -334,12 +368,16 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.CONTROL' include 'COMMON.SAXS' include 'COMMON.TORCNSTR' + include 'COMMON.TIME1' double precision evdw,evdw1,evdw2,evdw2_14,ees,eel_loc, & eello_turn3,eello_turn4,edfadis,estr,ehpb,ebe,ethetacnstr, & escloc,etors,edihcnstr,etors_d,esccor,ecorr,ecorr5,ecorr6,eturn6, & eliptran,Eafmforce,Etube, & esaxs_constr,ehomology_constr,edfator,edfanei,edfabet integer i,n_corr,n_corr1 +#ifdef TIMING_ENE + double precision time01 +#endif c write(iout,'(a,i2)')'Calling etotal_short ipot=',ipot c call flush(iout) if (modecalc.eq.12.or.modecalc.eq.14) then @@ -432,13 +470,13 @@ C FG slaves receive the WEIGHTS array wsaxs=weights(26) endif c write (iout,*),"Processor",myrank," BROADCAST weights" - call MPI_Bcast(c(1,1),maxres6,MPI_DOUBLE_PRECISION, + call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION, & king,FG_COMM,IERR) c write (iout,*) "Processor",myrank," BROADCAST c" - call MPI_Bcast(dc(1,1),maxres6,MPI_DOUBLE_PRECISION, + call MPI_Bcast(dc(1,1),6*nres,MPI_DOUBLE_PRECISION, & king,FG_COMM,IERR) c write (iout,*) "Processor",myrank," BROADCAST dc" - call MPI_Bcast(dc_norm(1,1),maxres6,MPI_DOUBLE_PRECISION, + call MPI_Bcast(dc_norm(1,1),6*nres,MPI_DOUBLE_PRECISION, & king,FG_COMM,IERR) c write (iout,*) "Processor",myrank," BROADCAST dc_norm" call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION, @@ -470,6 +508,9 @@ c call int_from_cart1(.false.) C C Compute the side-chain and electrostatic interaction energy C +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif goto (101,102,103,104,105,106) ipot C Lennard-Jones potential. 101 call elj_short(evdw) @@ -494,16 +535,31 @@ C C Calculate electrostatic (H-bonding) energy of the main chain. C 107 continue +#ifdef TIMING_ENE + time_evdw_short=time_evdw_short+MPI_Wtime()-time01 +#endif c c Calculate the short-range part of Evdwpp c +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif call evdwpp_short(evdw1) +#ifdef TIMING_ENE + time_eelec_short=time_eelec_short+MPI_Wtime()-time01 +#endif c c Calculate the short-range part of ESCp c +#ifdef TIMING_ENE + time01=MPI_Wtime() +#endif if (ipot.lt.6) then call escp_short(evdw2,evdw2_14) endif +#ifdef TIMING_ENE + time_escp_short=time_escp_short+MPI_Wtime()-time01 +#endif c c Calculate the bond-stretching energy c @@ -531,7 +587,13 @@ C energy function C C Calculate the SC local energy. C +#ifdef TIMING + time01=MPI_Wtime() +#endif call vec_and_deriv +#ifdef TIMING + time_vec=time_vec+MPI_Wtime()-time01 +#endif call esc(escloc) C C Calculate the virtual-bond torsional energy. diff --git a/source/unres/src-HCD-5D/gen_rand_conf.F b/source/unres/src-HCD-5D/gen_rand_conf.F index ea009b6..603e6c0 100644 --- a/source/unres/src-HCD-5D/gen_rand_conf.F +++ b/source/unres/src-HCD-5D/gen_rand_conf.F @@ -13,7 +13,7 @@ C Generate random conformation or chain cut and regrowth. logical overlap,back,fail integer nstart integer i,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit - double precision gen_theta,gen_phi,dist + double precision gen_theta,gen_phi,dist,ran_number cd print *,' CG Processor',me,' maxgen=',maxgen maxsi=100 cd write (iout,*) 'Gen_Rand_conf: nstart=',nstart @@ -60,12 +60,16 @@ c write(iout,*)'theta(3)=',rad2deg*theta(3) it1=iabs(itype(i-1)) it2=iabs(itype(i-2)) it=iabs(itype(i)) -c print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2, -c & ' nit=',nit,' niter=',niter,' maxgen=',maxgen + if (it.eq.ntyp1 .and. it1.eq.ntyp1) then + vbld(i)=ran_number(3.8d0,10.0d0) + vbld_inv(i)=1.0d0/vbld(i) + endif +c write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1, +c & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen phi(i+1)=gen_phi(i+1,it1,it) if (back) then phi(i)=gen_phi(i+1,it2,it1) -c print *,'phi(',i,')=',phi(i) +c write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it theta(i-1)=gen_theta(it2,phi(i-1),phi(i)) if (it2.ne.10 .and. it2.ne.ntyp1) then nsi=0 @@ -78,7 +82,12 @@ c print *,'phi(',i,')=',phi(i) endif call locate_next_res(i-1) endif - theta(i)=gen_theta(it1,phi(i),phi(i+1)) + if (it1.ne.ntyp1) then + theta(i)=gen_theta(it1,phi(i),phi(i+1)) + else + theta(i)=ran_number(1.326d0,2.548d0) + endif +c write (iout,*) "i",i," theta",theta(i) if (it1.ne.10 .and. it1.ne.ntyp1) then nsi=0 fail=.true. @@ -86,10 +95,12 @@ c print *,'phi(',i,')=',phi(i) call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail) nsi=nsi+1 enddo +c write (iout,*) "After forward SC generation:",nsi,maxsi if (nsi.gt.maxsi) return1 endif call locate_next_res(i) if (overlap(i-1)) then +c write (iout,*) "overlap",i-1 if (nit.lt.maxnit) then back=.true. nit=nit+1 @@ -107,6 +118,7 @@ c print *,'phi(',i,')=',phi(i) endif endif else +c write (iout,*) "No overlap",i-1 back=.false. nit=0 i=i+1 @@ -211,11 +223,55 @@ c-------------------------------------------------------------------------- double precision function gen_phi(i,it1,it2) implicit real*8 (a-h,o-z) include 'DIMENSIONS' + include 'COMMON.IOUNITS' include "COMMON.TORCNSTR" include 'COMMON.GEO' include 'COMMON.BOUNDS' - if (raw_psipred .or. ndih_constr.eq.0) then + include 'COMMON.INTERACT' + double precision sumprob(3) + double precision pinorm + external pinorm + if (ndih_constr.eq.0) then gen_phi=ran_number(-pi,pi) + else if (raw_psipred) then + if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 + & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then + ii=iconstr_dih(i) + sumprob(1)=vpsipred(2,ii) + sumprob(2)=sumprob(1)+vpsipred(3,ii) + sumprob(3)=sumprob(2)+vpsipred(1,ii) + aux=ran_number(0.0d0,sumprob(3)) +#ifdef DEBUG + write(iout,*)"gen_phi: residue i",i," ii",ii," vpsipred", + & vpsipred(2,ii),vpsipred(3,ii),vpsipred(1,ii)," sumprob", + & sumprob(1),sumprob(2),sumprob(3) + write (iout,*) "aux",aux +#endif + if (aux.le.sumprob(1)) then +#ifdef DEBUG + write (iout,*) "hel:", + & (phibound(1,i)-sdihed(1,ndih_constr))*rad2deg, + & (phibound(1,i)+sdihed(1,ndih_constr))*rad2deg +#endif + gen_phi=ran_number(phibound(1,i)-sdihed(1,ndih_constr) + & ,phibound(1,i)+sdihed(1,ndih_constr)) + else if (aux.le.sumprob(2)) then +#ifdef DEBUG + write (iout,*) "ext:", + & (phibound(2,i)-sdihed(2,ndih_constr))*rad2deg, + & (phibound(2,i)+sdihed(2,ndih_constr))*rad2deg +#endif + gen_phi=pinorm(ran_number(phibound(2,i)-sdihed(2,ndih_constr) + & ,phibound(2,i)+sdihed(2,ndih_constr))) + else +#ifdef DEBUG + write (iout,*) "coil:",-180.0,180.0 +#endif + gen_phi=ran_number(-pi,pi) + endif + else + gen_phi=ran_number(-pi,pi) + endif else C 8/13/98 Generate phi using pre-defined boundaries gen_phi=ran_number(phibound(1,i),phibound(2,i)) @@ -288,7 +344,7 @@ c------------------------------------------------------------------------- fail=.false. if (the.eq.0.0D0 .or. the.eq.pi) then #ifdef MPI - write (*,'(a,i4,a,i3,a,1pe14.5)') + write (iout,'(a,i4,a,i3,a,1pe14.5)') & 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the #else cd write (iout,'(a,i3,a,1pe14.5)') @@ -778,23 +834,24 @@ c overlapping residues left, or false otherwise (success) include 'COMMON.IOUNITS' logical had_overlaps,fail,scfail integer ioverlap(maxres),ioverlap_last + integer maxit_corr /5000/ had_overlaps=.false. - call overlap_sc_list(ioverlap,ioverlap_last) + call overlap_sc_list(ioverlap,ioverlap_last,.false.) if (ioverlap_last.gt.0) then write (iout,*) '#OVERLAPing residues ',ioverlap_last - write (iout,'(18i5)') (ioverlap(k),k=1,ioverlap_last) + write (iout,'(15i6)') (ioverlap(k),k=1,ioverlap_last) had_overlaps=.true. endif maxsi=1000 - do k=1,1000 + do k=1,maxit_corr if (ioverlap_last.eq.0) exit do ires=1,ioverlap_last i=ioverlap(ires) iti=iabs(itype(i)) - if (iti.ne.10) then + if (iti.ne.10 .and. iti.lt.ntyp1) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) @@ -805,16 +862,15 @@ c overlapping residues left, or false otherwise (success) if(fail) goto 999 endif enddo - c write (iout,*) "before chaincuild overlap_sc_list: dc0",dc(:,0) c call chainbuild_extconf c write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0) - call overlap_sc_list(ioverlap,ioverlap_last) - write (iout,*) 'Overlaping residues ',ioverlap_last, - & (ioverlap(j),j=1,ioverlap_last) + call overlap_sc_list(ioverlap,ioverlap_last,.false.) + write (iout,*)'#Overlaping residues @iter',k,":",ioverlap_last + write (iout,*)'Residue list:',(ioverlap(j),j=1,ioverlap_last) enddo - if (k.le.1000.and.ioverlap_last.eq.0) then + if (k.le.maxit_corr.and.ioverlap_last.eq.0) then scfail=.false. if (had_overlaps) then write (iout,*) '#OVERLAPing all corrected after ',k, @@ -823,22 +879,29 @@ c write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0) else scfail=.true. write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last - write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last) + write (iout,'(15i6)') (ioverlap(j),j=1,ioverlap_last) endif return 999 continue - write (iout,'(a30,i5,a12,i4)') + write (iout,'(a30,i5,a12,i6)') & '#OVERLAP FAIL in gen_side after',maxsi, & 'iter for RES',i scfail=.true. return end - subroutine overlap_sc_list(ioverlap,ioverlap_last) - implicit real*8 (a-h,o-z) + subroutine overlap_sc_list(ioverlap,ioverlap_last,lprn) + implicit none include 'DIMENSIONS' +#ifdef MPI + include "mpif.h" + include "COMMON.SETUP" + integer ierror + integer ioverlap_last_tab(0:max_fg_procs-1), + & ioverlap_all(maxres*max_fg_procs),displs(0:max_fg_procs-1) +#endif include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.IOUNITS' @@ -847,19 +910,47 @@ c write (iout,*) "after chaincuild overlap_sc_list: dc0",dc(:,0) include 'COMMON.FFIELD' include 'COMMON.VAR' include 'COMMON.CALC' - logical fail + integer ii,itypi,itypj,itypi1,ind,ikont + logical fail,lprn integer ioverlap(maxres),ioverlap_last - data redfac /0.5D0/ + double precision redfac /0.5D0/ + double precision rrij,rij_shift,sig0ij,xi,yi,zi,rcomp,sig + double precision dist - write (iout,*) "overlap_sc_list" +#ifdef MPI + if (nfgtasks.gt.1) then + if (fg_rank.eq.0) + & call MPI_Bcast(11,1,MPI_INTEGER,king,FG_COMM,IERROR) + call MPI_Bcast(c(1,1),6*nres,MPI_DOUBLE_PRECISION,king,FG_COMM, + & IERROR) + call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,king, + & FG_COMM,IERROR) + call MPI_Bcast(dc_norm(1,0),6*(nres+1),MPI_DOUBLE_PRECISION, + & king,FG_COMM,IERROR) + endif +#endif +c write (iout,*) "overlap_sc_list" c write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e - write(iout,*) "nnt",nnt," nct",nct +c write(iout,*) "nnt",nnt," nct",nct ioverlap_last=0 C Check for SC-SC overlaps and mark residues c print *,'>>overlap_sc nnt=',nnt,' nct=',nct ind=0 +#ifdef DEBUG + write (iout,*) "FG proecssor",fg_rank," g_listscsc_start", + & g_listscsc_start," g_listscsc_end",g_listscsc_end + write (*,*) "FG proecssor",fg_rank," g_listscsc_start", + & g_listscsc_start," g_listscsc_end",g_listscsc_end +#endif c do i=iatsc_s,iatsc_e - do i=nnt,nct + do ikont=g_listscsc_start,g_listscsc_end +c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last", +c & ioverlap_last + i=newcontlisti(ikont) + j=newcontlistj(ikont) +c do i=nnt,nct +c write (*,*) "FG processor",fg_rank," loop begins ioverlap_last", +c & ioverlap_last,"ikont i j",ikont,i,j itypi=iabs(itype(i)) itypi1=iabs(itype(i+1)) if (itypi.eq.ntyp1) cycle @@ -873,7 +964,7 @@ c do i=iatsc_s,iatsc_e c c do iint=1,nint_gr(i) c do j=istart(i,iint),iend(i,iint) - do j=i+1,nct +c do j=i+1,nct ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -919,10 +1010,16 @@ c & " sig",sig," sig0ij",sig0ij c if ( rij_shift.le.0.0D0 ) then if ( rij_shift/sig0ij.le.0.1D0 ) then c write (iout,*) "overlap",i,j - write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)') - & 'overlap SC-SC: i=',i,' j=',j, - & ' dist=',dist(nres+i,nres+j),' rcomp=', - & rcomp,1.0/rij,rij_shift + if (lprn) then + write (iout,'(a,i5,a,i5,a,f10.5,a,3f10.5)') + & 'overlap SC-SC: i=',i,' j=',j, + & ' dist=',dist(nres+i,nres+j),' rcomp=', + & rcomp,1.0/rij,rij_shift + write (*,'(a,i2,a,i5,a,i5,a,f10.5,a,3f10.5)') + & 'FG processor',fg_rank,' overlap SC-SC: i=',i,' j=',j, + & ' dist=',dist(nres+i,nres+j),' rcomp=', + & rcomp,1.0/rij,rij_shift + endif ioverlap_last=ioverlap_last+1 ioverlap(ioverlap_last)=i do k=1,ioverlap_last-1 @@ -933,9 +1030,73 @@ c write (iout,*) "overlap",i,j do k=1,ioverlap_last-1 if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1 enddo +c write(*,*) "FG processor",fg_rank,i,j," ioverlap_last", +c & ioverlap_last," ioverlap",(ioverlap(k),k=1,ioverlap_last) endif - enddo +c enddo c enddo enddo +#ifdef MPI +#ifdef DEBUG + write (iout,*) "FG Processor",fg_rank," ioverlap_last", + & ioverlap_last," ioverlap",(ioverlap(i),i=1,ioverlap_last) + write (*,*) "FG Processor",fg_rank," ioverlap_last",ioverlap_last, + & " ioverlap",(ioverlap(i),i=1,ioverlap_last) + call flush(iout) +#endif + if (nfgtasks.eq.1) return +#ifdef DEBUG + write (iout,*) "Before MPI_Gather" + call flush(iout) +#endif + call MPI_Gather(ioverlap_last,1,MPI_INTEGER,ioverlap_last_tab, + & 1,MPI_INTEGER,king,FG_COMM,IERROR) +#ifdef DEBUG + write (iout,*) "After MPI_Gather" + call flush(iout) +#endif +#ifdef DEBUG + if (myrank.eq.king) + & write (iout,*) "FG Processor",fg_rank,"ioverlap_last_tab", + & (ioverlap_last_tab(i),i=0,nfgtasks-1) + call flush(iout) +#endif + displs(0)=0 + do i=1,nfgtasks-1 + displs(i)=displs(i-1)+ioverlap_last_tab(i-1) + enddo + call MPI_Gatherv(ioverlap,ioverlap_last,MPI_INTEGER, + & ioverlap_all,ioverlap_last_tab,displs,MPI_INTEGER,king, + & FG_COMM,IERROR) +#ifdef DEBUG + write (iout,*) "After Gatherv" + call flush(iout) +#endif + if (fg_rank.gt.0) return + ioverlap_last=0 + do i=0,nfgtasks-1 + ioverlap_last=ioverlap_last+ioverlap_last_tab(i) + enddo +#ifdef DEBUG + write (iout,*) "ioverlap_last",ioverlap_last," ioverlap_last", + & (ioverlap_all(i),i=1,ioverlap_last) + call flush(iout) +#endif + ii=0 + do i=1,ioverlap_last + ioverlap(ii+1)=ioverlap_all(i) + do j=ii,1,-1 + if (ioverlap(ii+1).eq.ioverlap(j)) goto 11 + enddo + ii=ii+1 + 11 continue + enddo + ioverlap_last=ii +#ifdef DEBUG + write (iout,*) "After summing: ioverlap_last",ioverlap_last, + & " ioverlap",(ioverlap(i),i=1,ioverlap_last) + call flush(iout) +#endif +#endif return end diff --git a/source/unres/src-HCD-5D/gen_rand_conf_mchain.F b/source/unres/src-HCD-5D/gen_rand_conf_mchain.F new file mode 100644 index 0000000..4614e87 --- /dev/null +++ b/source/unres/src-HCD-5D/gen_rand_conf_mchain.F @@ -0,0 +1,424 @@ + subroutine gen_rand_conf_mchain(nstart0,*) +C Generate random conformation or chain cut and regrowth. + implicit none + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.LOCAL' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.MCM' + include 'COMMON.GEO' + include 'COMMON.CONTROL' + logical overlap_mchain,back,fail + integer nstart,nstart0 + integer i,ii,iii,j,k,it,it1,it2,nit,niter,nsi,maxsi,maxnit + integer igr,iequi,ichain,nnres + double precision aux + double precision gen_theta,gen_phi,dist,ran_number,scalar +c write (iout,*) 'gen_rand_conf_mchain: maxgen=',maxgen + nstart=nstart0 + maxsi=100 +c write (iout,*) 'Gen_Rand_conf_mchain: nstart=',nstart, +c & " nchain_group",nchain_group + + DO IGR=1,NCHAIN_GROUP + + DO IEQUI=1,NEQUIV(IGR) + + ichain=iequiv(iequi,igr) + + i=chain_border1(1,ichain)+nstart-1 + if (nstart.eq.1) then + do j=1,3 + c(j,i)=ran_number(-15.0d0,15.0d0) + dc(j,i-1)=c(j,i) + enddo + endif + if (nstart.le.2) then + + do j=1,3 + dc_norm(j,i)=ran_number(-1.0d0,1.0d0) + enddo + aux=dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))) + do j=1,3 + dc_norm(j,i)=dc_norm(j,i)/aux + enddo + if (itype(i).eq.ntyp1) then + do j=1,3 + dc(j,i)=1.9d0*dc_norm(j,i) + enddo + else + do j=1,3 + dc(j,i)=3.8d0*dc_norm(j,i) + enddo + endif + do j=1,3 + c(j,i+1)=c(j,i)+dc(j,i) + enddo + endif +c if (nstart.lt.5) then + if (nstart.le.2) then + it1=iabs(itype(i+2)) + phi(i+3)=gen_phi(i+3,iabs(itype(i+1)),iabs(itype(i+2))) +c write(iout,*)'phi(4)=',rad2deg*phi(4) + theta(i+2)=gen_theta(iabs(itype(i+2)),pi,phi(i+3)) + if (it1.ne.10) then + nsi=0 + fail=.true. + do while (fail.and.nsi.le.maxsi) + call gen_side(it1,theta(i+2),alph(i+1),omeg(i+1),fail) + nsi=nsi+1 + enddo + if (nsi.gt.maxsi) then + write (iout,'(a,i7,a,i7,a,i7,a,i7)') + & 'Problem in SC rotamer generation, residue',ii, + & ' chain',ichain,' chain group',igr,'. Increase MAXSI.' + return1 + endif + endif ! it1.ne.10 + call orig_frame_chain(i) + else +c call orig_frame_chain(i-1) +c write (iout,*) "calling refsys",i + call refsys(i,i-1,i-2,prod(1,1,i-1), + & prod(1,2,i-1),prod(1,3,i-1),fail) +c write (iout,*) "after refsys",i +#ifdef DEBUG + write (iout,*) "dc_norm(:",i-1,") and prod" + do j=1,3 + write (iout,*) j,dc_norm(j,i-1),(prod(j,k,i-1),k=1,3) + enddo +#endif + endif + + ENDDO + + nstart=nstart+1 + ii=nstart + + maxnit=5000 + + nit=0 + niter=0 + back=.false. + nnres=chain_border1(2,iequiv(1,igr))- + & chain_border1(1,iequiv(1,igr))+1 +#ifdef DEBUG + write (iout,*) "chain group",igr," chains", + & (iequiv(j,igr),j=1,nequiv(igr)) + write (iout,*) "ii",ii," nnres",nnres +#endif + do while (ii.le. nnres .and. niter.lt.maxgen) + + ichain=iequiv(1,igr) + i=ii-1+chain_border1(1,ichain) +#ifdef DEBUG + write (iout,*) "ii",ii," nnres",nnres," ichain",ichain," i",i, + & "niter",niter," back",back," nstart",nstart +#endif + if (ii.lt.nstart) then +c if(iprint.gt.1) then + write (iout,'(/80(1h*)/2a/80(1h*))') + & 'Generation procedure went down to ', + & 'chain beginning. Cannot continue...' + write (*,'(/80(1h*)/2a/80(1h*))') + & 'Generation procedure went down to ', + & 'chain beginning. Cannot continue...' +c endif + return1 + endif + it1=iabs(itype(i-1)) + it2=iabs(itype(i-2)) + it=iabs(itype(i)) + if (it.eq.ntyp1 .and. it1.eq.ntyp1) then + vbld(i)=ran_number(3.8d0,10.0d0) + vbld_inv(i)=1.0d0/vbld(i) + endif +#ifdef DEBUG + write (iout,*) 'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1, + & ' it2=',it2,' nit=',nit,' niter=',niter,' maxgen=',maxgen +#endif + phi(i)=gen_phi(i,it2,it1) + phi(i+1)=gen_phi(i+1,it1,it) +#ifdef DEBUG + write (iout,*) "phi",i,phi(i)," phi",i+1,phi(i+1) +#endif + do iequi=2,nequiv(igr) + iii=ii+chain_border1(1,iequiv(iequi,igr)) + phi(iii)=phi(i+1) + enddo +#ifdef CHUJ + if (back) then + phi(i)=gen_phi(i+1,it2,it1) +#ifdef DEBUG + write(iout,*) 'phi(',i,')=',phi(i)," type",it1,it2,it +#endif + theta(i-1)=gen_theta(it2,phi(i-1),phi(i)) + if (theta(i-1).gt.2.68780478d0) theta(i-1)=2.68780478d0 + do iequi=2,nequiv(igr) + iii=ii-1+chain_border1(1,iequiv(iequi,igr)) + phi(iii)=phi(i+1) + theta(iii-1)=theta(i-1) + enddo + if (it2.ne.10 .and. it2.ne.ntyp1) then + nsi=0 + fail=.true. + do while (fail.and.nsi.le.maxsi) + call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail) + nsi=nsi+1 + enddo + do iequi=2,nequiv(igr) + iii=ii-3+chain_border1(1,iequiv(iequi,igr)) + alph(iii)=alph(i-2) + omeg(iii)=omeg(i-2) + enddo +#ifdef DEBUG + write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail +#endif + if (nsi.gt.maxsi) then + write (iout,'(a,i7,a,i7,a,i7,a,i7)') + & 'Problem in SC rotamer generation, residue',ii, + & ' chain',ichain,' chain group',igr,'. Increase MAXSI.' + return1 + endif + + endif +c call locate_next_res(i-1) + endif +#endif + if (it1.ne.ntyp1) then + theta(i)=gen_theta(it1,phi(i),phi(i+1)) + if (theta(i).gt.2.68780478d0) theta(i)=2.68780478d0 + else + theta(i)=ran_number(1.326d0,2.548d0) + endif +#ifdef DEBUG + write (iout,*) "ii",ii," i",i," it1",it1, + & " theta",theta(i)," phi", + & phi(i),phi(i+1) +#endif + if (it1.ne.10 .and. it1.ne.ntyp1) then + nsi=0 + fail=.true. + do while (fail.and.nsi.le.maxsi) +c theta(i)=gen_theta(it1,phi(i),phi(i+1)) + call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail) + nsi=nsi+1 + enddo +#ifdef DEBUG + write (iout,*) "alpha",alph(i)," omeg",omeg(i)," fail",fail +#endif +c write (iout,*) "After forward SC generation:",nsi,maxsi + if (nsi.gt.maxsi) then + write (iout,'(a,i7,a,i7,a,i7,a,i7)') + & 'Problem in SC rotamer generation, residue',ii, + & ' chain',ichain,' chain group',igr,'. Increase MAXSI.' + return1 + endif + + endif + + do iequi=2,nequiv(igr) + iii=ii-1+chain_border1(1,iequiv(iequi,igr)) + theta(iii)=theta(i) + phi(iii)=phi(i) + alph(iii-1)=alph(i-1) + omeg(iii-1)=omeg(i-1) + enddo + + DO IEQUI=1,NEQUIV(IGR) + + ichain=iequiv(iequi,igr) + i=ii-1+chain_border1(1,ichain) +#ifdef CHUJ + if (back) call locate_next_res(i-1) +#endif + call locate_next_res(i) +#ifdef DEBUG + write (iout,*) "i",i," ii",ii," ichain",ichain + write (iout,*) theta(i)*rad2deg,phi(i)*rad2deg, + & alph(i-1)*rad2deg,omeg(i-1)*rad2deg + write (iout,*) (c(j,i),j=1,3),(c(j,i+nres-1),j=1,3) +#endif + if (overlap_mchain(i-1,ii-1,ichain,igr)) then +#ifdef DEBUG + write (iout,*) "***********overlap",i-1," nit",nit +#endif + if (nit.lt.maxnit) then + back=.true. + nit=nit+1 + exit + else +#ifdef DEBUG + write (iout,*) "***********overlap maxnit exceeded",nit +#endif + nit=0 + if (ii.gt.3) then + back=.true. + ii=ii-1 +c write (iout,*) "ii",ii + exit + else + write (iout,'(a,i7,a,i7,a,i7,a,i7)') + & 'Cannot generate non-overlaping conformation, residue',ii, + & ' chain',ichain,' chain group',igr,'. Increase MAXNIT.' + return1 + endif + endif + else +c write (iout,*) "No overlap",i-1 + back=.false. +c nit=0 +c i=i+1 + endif + + ENDDO + + if (.not.back) then +#ifdef DEBUG + write (iout,*) "++++++++++Successful generation",igr,ichain,ii +#endif + ii=ii+1 + nit=0 + endif + back=.false. +c write (iout,*) "ii",ii + niter=niter+1 + + enddo + if (niter.ge.maxgen) then + write (iout,'(a,2i7,a,i7,a,i7,a,i7)') + & 'Too many trials in conformation generation',niter,maxgen, + & ' chain group',igr,' chain',ichain,' residue',ii + return1 + endif + + ENDDO + + do i=2,nres + if (itype(i).eq.ntyp1) then + do j=1,3 + dc(j,i-1)=c(j,i)-c(j,i-1) + enddo + endif + enddo + + return + end +c------------------------------------------------------------------------- + logical function overlap_mchain(i,ii,ichain,igr) + implicit none + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + double precision redfac /0.5D0/,redfacp /0.8d0/,redfacscp /0.8d0/ + integer i,ii,ichain,igr,j,jj,jchain,jgr,k,iti,itj,iteli,itelj + double precision rcomp + double precision dist + logical lprn /.false./ + overlap_mchain=.false. + iti=iabs(itype(i)) + if (iti.gt.ntyp) return +C Check for SC-SC overlaps. +cd print *,'nnt=',nnt,' nct=',nct +#ifdef DEBUG + write (iout,*) "overlap_mchain i",i," ii",ii," ichain",ichain, + & " igr",igr," iti",iti +#endif + do j=nnt,nct + if (itype(j).eq.ntyp1) cycle + jchain=ireschain(j) +c write (iout,*) "overlap_mchain j",j," jj",jj," jchain",jchain, +c & " itype",itype(j) + jgr=mapchain(jchain) + jj=j-chain_border1(1,jchain)+1 + itj=iabs(itype(j)) +c write (iout,*) "jgr",jgr + if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-1 + & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle + if (j.lt.i-1 .or. ipot.ne.4) then + rcomp=sigmaii(iti,itj) + else + rcomp=sigma(iti,itj) + endif +cd print *,'j=',j + if (dist(nres+i,nres+j).lt.redfac*rcomp) then + overlap_mchain=.true. + if (lprn) write(iout,*)'overlap_mchain, SC-SC: i=',i,' j=',j, + & ' ichain',ichain,' jchain',jchain, + & ' dist=',dist(nres+i,nres+j),' rcomp=', + & rcomp*redfac + return + endif + enddo +#ifdef CHUJ +C Check for overlaps between the added peptide group and the preceding +C SCs. + iteli=itel(i) + if (iteli.gt.0) then + do j=1,3 + c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1)) + enddo + do j=nnt,nct + itj=iabs(itype(j)) + if (itj.eq.ntyp1) cycle + jchain=ireschain(j) + jgr=mapchain(jchain) + jj=j-chain_border1(1,jchain)+1 + if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2 + & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle +cd print *,'overlap_mchain, p-Sc: i=',i,' j=',j, +cd & ' dist=',dist(nres+j,maxres2+1) + if (dist(j,maxres2+1).lt.4.0D0*redfacscp) then + if (lprn) write (iout,*) 'overlap_mchain, p-SC: i=',i,' j=',j, + & ' ichain',ichain,' jchain',jchain, + & ' dist=',dist(nres+j,maxres2+1),' rcomp=', + & 4.0d0*redfac + overlap_mchain=.true. + return + endif + enddo + endif +C Check for p-p overlaps +#endif + iteli=itel(i) + if (iteli.eq.0) return + do j=1,3 + c(j,maxres2+1)=0.5D0*(c(j,i)+c(j,i+1)) + enddo +c do j=nnt,i-2 + do j=nnt,nct + itelj=itel(j) + if (itelj.eq.0) cycle + do k=1,3 + c(k,maxres2+2)=0.5D0*(c(k,j)+c(k,j+1)) + enddo + jchain=ireschain(j) + jgr=mapchain(jchain) +c if(iteli.ne.0.and.itelj.ne.0)then +c write (iout,*) i,j,dist(maxres2+1,maxres2+2),rpp(iteli,itelj) + jj=j-chain_border1(1,jchain)+1 +c write (iout,*) "jgr",jgr + if(igr.eq.jgr.and.jj.gt.ii .or. ichain.eq.jchain .and. j.gt.i-2 + & .or. jgr.gt.igr .and. jj.gt.nran_start) cycle +#ifdef DEBUG + write (iout,*)'overlap_mchain, p-p: i=',i,' j=',j,' igr',igr, + & ' jgr',jgr,' ichain',ichain,' jchain',jchain, + & ' dist=',dist(maxres2+1,maxres2+2) +#endif + if (dist(maxres2+1,maxres2+2).lt.rpp(iteli,itelj)*redfacp) then + if (lprn) write (iout,*) 'overlap_mchain, p-p: i=',i,' j=',j, + & ' ichain',ichain,' jchain',jchain, + & ' dist=',dist(maxres2+1,maxres2+2),' rcomp=', + & rpp(iteli,itelj)*redfacp + overlap_mchain=.true. + return + endif +c endif + enddo + return + end + diff --git a/source/unres/src-HCD-5D/geomout.F b/source/unres/src-HCD-5D/geomout.F index d1a3a87..3dcde10 100644 --- a/source/unres/src-HCD-5D/geomout.F +++ b/source/unres/src-HCD-5D/geomout.F @@ -16,7 +16,10 @@ #endif character*50 tytul integer iunit - character*1 chainid(10) /'A','B','C','D','E','F','G','H','I','J'/ + character*1 chainid(52) /'A','B','C','D','E','F','G','H','I','J', + & 'K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z', + & 'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p', + & 'q','r','s','t','u','v','w','x','y','z'/ integer ica(maxres) integer i,j,k,iti,itj,itk,itl,iatom,ichain,ires double precision etot @@ -90,8 +93,8 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 do i=1,nss if (dyn_ss) then write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') - & 'SSBOND',i,'CYS',idssb(i)-nnt+1, - & 'CYS',jdssb(i)-nnt+1 + & 'SSBOND',i,'CYS',iss(idssb(i))-nnt+1, + & 'CYS',iss(jdssb(i))-nnt+1 else write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') & 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres, @@ -107,9 +110,10 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 iti=itype(i) if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then ichain=ichain+1 + if (ichain.gt.52) ichain=1 ires=0 write (iunit,'(a)') 'TER' - else + else if (iti.ne.ntyp1) then ires=ires+1 iatom=iatom+1 ica(i)=iatom @@ -149,7 +153,7 @@ cmodel write (iunit,'(a5,i6)') 'MODEL',1 write (iunit,'(a6)') 'ENDMDL' 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3) 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3) - 30 FORMAT ('CONECT',8I5) + 30 FORMAT ('CONECT',8I7) return end c------------------------------------------------------------------------------ @@ -541,7 +545,7 @@ C print *,'A CHUJ',potEcomp(23) if(itime.eq.0) then write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2, & ",31a12)" - write (istat,format) "#","", + write (istat,format) "#"," ", & (ename(print_order(i)),i=1,nprint_ene) endif write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2, diff --git a/source/unres/src-HCD-5D/gradient_p.F b/source/unres/src-HCD-5D/gradient_p.F index 67275ed..af7978b 100644 --- a/source/unres/src-HCD-5D/gradient_p.F +++ b/source/unres/src-HCD-5D/gradient_p.F @@ -210,6 +210,7 @@ C------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.TIME1' integer i,j,kk + double precision time00,time01 c c This subrouting calculates total Cartesian coordinate gradient. c The subroutine chainbuild_cart and energy MUST be called beforehand. diff --git a/source/unres/src-HCD-5D/initialize_p.F b/source/unres/src-HCD-5D/initialize_p.F index 710f907..6a297b8 100644 --- a/source/unres/src-HCD-5D/initialize_p.F +++ b/source/unres/src-HCD-5D/initialize_p.F @@ -351,7 +351,7 @@ c------------------------------------------------------------------------- 1 "WSC ","WSCP ","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC", ! 8 9 10 11 12 13 14 8 "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR ","WTORD ", -!  15 16 17 18 19 20 21 +!  15 16 17 18 19 20 21 5 "WSTRAIN","WVDWPP","WBOND","SCAL14","WDIHC","WUMB","WSCCOR", ! 22 23 24 25 26 27 28 2 "WLT","WAFM","WTHETCNSR","WTUBE","WSAXS","WHOMO","WDFAD", @@ -453,8 +453,13 @@ c--------------------------------------------------------------------------- integer iturn3_start_all,iturn3_end_all,iturn4_start_all, & iturn4_end_all,iatel_s_all, & iatel_e_all,ielstart_all,ielend_all,ntask_cont_from_all, - & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all, - & n_sc_int_tot,my_sc_inds,my_sc_inde,ind_sctint,ind_scint_old + & itask_cont_from_all,ntask_cont_to_all,itask_cont_to_all + integer*8 n_sc_int_tot,my_sc_inds,my_sc_inde,ind_scint, + & ind_scint_old,nele_int_tot,ind_eleint,my_ele_inds,my_ele_inde, + & ind_eleint_old,nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw, + & ind_eleint_vdw,ind_eleint_vdw_old,nscp_int_tot,my_scp_inds, + & my_scp_inde,ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start, + & ngrad_end common /przechowalnia/ iturn3_start_all(0:max_fg_procs), & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs), & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs), @@ -466,12 +471,8 @@ c--------------------------------------------------------------------------- & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1) integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP logical scheck,lprint,flag - integer i,j,k,ii,jj,iint,npept,nele_int_tot,ind_eleint,ind_scint, - & my_ele_inds,my_ele_inde,ind_eleint_old,nele_int_tot_vdw, - & my_ele_inds_vdw,my_ele_inde_vdw,ind_eleint_vdw,ijunk, - & ind_eleint_vdw_old,nscp_int_tot,my_scp_inds,my_scp_inde, - & ind_scpint,ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end, - & iaux,ind_typ,ncheck_from,ncheck_to,ichunk + integer i,j,k,ii,jj,iint,npept, + & ijunk,iaux,ind_typ,ncheck_from,ncheck_to,ichunk #ifdef MPI integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs), & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs) @@ -486,14 +487,14 @@ C... to deal with by current processor. lprint=energy_dec if (lprint) &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct - n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss - call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde) + n_sc_int_tot=int(nct-nnt+1,8)*int(nct-nnt,8)/2-nss + call int_bounds8(n_sc_int_tot,my_sc_inds,my_sc_inde) if (lprint) & write (iout,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',MyRank, & ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds, & ' my_sc_inde',my_sc_inde - ind_sctint=0 + ind_scint=0 iatsc_s=0 iatsc_e=0 #endif @@ -530,7 +531,7 @@ cd write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj if (jj.eq.i+1) then #ifdef MPI c write (iout,*) 'jj=i+1' - call int_partition(ind_scint,my_sc_inds,my_sc_inde,i, + call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i, & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12) #else nint_gr(i)=1 @@ -540,7 +541,7 @@ c write (iout,*) 'jj=i+1' else if (jj.eq.nct) then #ifdef MPI c write (iout,*) 'jj=nct' - call int_partition(ind_scint,my_sc_inds,my_sc_inde,i, + call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i, & iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12) #else nint_gr(i)=1 @@ -549,10 +550,10 @@ c write (iout,*) 'jj=nct' #endif else #ifdef MPI - call int_partition(ind_scint,my_sc_inds,my_sc_inde,i, + call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i, & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12) ii=nint_gr(i)+1 - call int_partition(ind_scint,my_sc_inds,my_sc_inde,i, + call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i, & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12) #else nint_gr(i)=2 @@ -564,7 +565,7 @@ c write (iout,*) 'jj=nct' endif else #ifdef MPI - call int_partition(ind_scint,my_sc_inds,my_sc_inde,i, + call int_partition8(ind_scint,my_sc_inds,my_sc_inde,i, & iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12) #else nint_gr(i)=1 @@ -588,9 +589,10 @@ c write (iout,*) 'jj=nct' & ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e #endif if (lprint) then + write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e write (iout,'(a)') 'Interaction array:' do i=iatsc_s,iatsc_e - write (iout,'(i3,2(2x,2i3))') + write (iout,'(i7,2(2x,2i7))') & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i)) enddo endif @@ -598,8 +600,8 @@ c write (iout,*) 'jj=nct' #ifdef MPI C Now partition the electrostatic-interaction array npept=nct-nnt - nele_int_tot=(npept-ispp)*(npept-ispp+1)/2 - call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde) + nele_int_tot=int(npept-ispp,8)*int(npept-ispp+1,8)/2 + call int_bounds8(nele_int_tot,my_ele_inds,my_ele_inde) if (lprint) & write (*,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',MyRank, @@ -611,14 +613,14 @@ C Now partition the electrostatic-interaction array ind_eleint_old=0 do i=nnt,nct-3 ijunk=0 - call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i, + call int_partition8(ind_eleint,my_ele_inds,my_ele_inde,i, & iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13) enddo ! i 13 continue if (iatel_s.eq.0) iatel_s=1 - nele_int_tot_vdw=(npept-2)*(npept-2+1)/2 + nele_int_tot_vdw=int(npept-2,8)*int(npept-2+1,8)/2 c write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw - call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw) + call int_bounds8(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw) c write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw, c & " my_ele_inde_vdw",my_ele_inde_vdw ind_eleint_vdw=0 @@ -627,7 +629,7 @@ c & " my_ele_inde_vdw",my_ele_inde_vdw iatel_e_vdw=0 do i=nnt,nct-3 ijunk=0 - call int_partition(ind_eleint_vdw,my_ele_inds_vdw, + call int_partition8(ind_eleint_vdw,my_ele_inds_vdw, & my_ele_inde_vdw,i, & iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i), & ielend_vdw(i),*15) @@ -655,15 +657,15 @@ c & " ielend_vdw",ielend_vdw(i) & ' absolute rank',MyRank write (iout,*) 'Electrostatic interaction array:' do i=iatel_s,iatel_e - write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i) + write (iout,'(i7,2(2x,2i7))') i,ielstart(i),ielend(i) enddo endif ! lprint c iscp=3 iscp=2 C Partition the SC-p interaction array #ifdef MPI - nscp_int_tot=(npept-iscp+1)*(npept-iscp+1) - call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde) + nscp_int_tot=int(npept-iscp+1,8)*int(npept-iscp+1,8) + call int_bounds8(nscp_int_tot,my_scp_inds,my_scp_inde) if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor, & ' absolute rank',myrank, & ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds, @@ -675,20 +677,20 @@ C Partition the SC-p interaction array do i=nnt,nct-1 if (i.lt.nnt+iscp) then cd write (iout,*) 'i.le.nnt+iscp' - call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i, + call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i, & iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1), & iscpend(i,1),*14) else if (i.gt.nct-iscp) then cd write (iout,*) 'i.gt.nct-iscp' - call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i, + call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i, & iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1), & iscpend(i,1),*14) else - call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i, + call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i, & iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1), & iscpend(i,1),*14) ii=nscp_gr(i)+1 - call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i, + call int_partition8(ind_scpint,my_scp_inds,my_scp_inde,i, & iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii), & iscpend(i,ii),*14) endif @@ -719,7 +721,7 @@ cd write (iout,*) 'i.gt.nct-iscp' if (lprint) then write (iout,'(a)') 'SC-p interaction array:' do i=iatscp_s,iatscp_e - write (iout,'(i3,2(2x,2i3))') + write (iout,'(i7,2(2x,2i7))') & i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) enddo endif ! lprint @@ -782,24 +784,24 @@ c & " ivec_start",ivec_start," ivec_end",ivec_end endif c nsumgrad=(nres-nnt)*(nres-nnt+1)/2 c nlen=nres-nnt+1 - nsumgrad=(nres-nnt)*(nres-nnt+1)/2 - nlen=nres-nnt+1 - call int_bounds(nsumgrad,ngrad_start,ngrad_end) - igrad_start=((2*nlen+1) - & -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2 - jgrad_start(igrad_start)= - & ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2 - & +igrad_start - jgrad_end(igrad_start)=nres - igrad_end=((2*nlen+1) - & -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2 - if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1 - jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2 - & +igrad_end - do i=igrad_start+1,igrad_end-1 - jgrad_start(i)=i+1 - jgrad_end(i)=nres - enddo +c nsumgrad=(nres-nnt)*(nres-nnt+1)/2 +c nlen=nres-nnt+1 +c call int_bounds(nsumgrad,ngrad_start,ngrad_end) +c igrad_start=((2*nlen+1) +c & -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2 +c jgrad_start(igrad_start)= +c & ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2 +c & +igrad_start +c jgrad_end(igrad_start)=nres +c igrad_end=((2*nlen+1) +c & -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2 +c if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1 +c jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2 +c & +igrad_end +c do i=igrad_start+1,igrad_end-1 +c jgrad_start(i)=i+1 +c jgrad_end(i)=nres +c enddo if (lprint) then write (*,*) 'Processor:',fg_rank,' CG group',kolor, & ' absolute rank',myrank, @@ -818,13 +820,13 @@ c nlen=nres-nnt+1 & ' ithetaconstr_start',ithetaconstr_start,' ithetaconstr_end', & ithetaconstr_end - write (*,*) 'Processor:',fg_rank,myrank,' igrad_start', - & igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start, - & ' ngrad_end',ngrad_end - do i=igrad_start,igrad_end - write(*,*) 'Processor:',fg_rank,myrank,i, - & jgrad_start(i),jgrad_end(i) - enddo +c write (*,*) 'Processor:',fg_rank,myrank,' igrad_start', +c & igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start, +c & ' ngrad_end',ngrad_end +c do i=igrad_start,igrad_end +c write(*,*) 'Processor:',fg_rank,myrank,i, +c & jgrad_start(i),jgrad_end(i) +c enddo endif if (nfgtasks.gt.1) then call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1, @@ -1502,6 +1504,31 @@ c--------------------------------------------------------------------------- return end c--------------------------------------------------------------------------- + subroutine int_bounds8(total_ints,lower_bound,upper_bound) + implicit none + include 'DIMENSIONS' + include 'mpif.h' + include 'COMMON.SETUP' + integer*8 total_ints,lower_bound,upper_bound,nint + integer*8 int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs) + integer i,nexcess + nint=total_ints/nfgtasks + do i=1,nfgtasks + int4proc(i-1)=nint + enddo + nexcess=total_ints-nint*nfgtasks + do i=1,nexcess + int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1 + enddo + lower_bound=0 + do i=0,fg_rank-1 + lower_bound=lower_bound+int4proc(i) + enddo + upper_bound=lower_bound+int4proc(fg_rank) + lower_bound=lower_bound+1 + return + end +c--------------------------------------------------------------------------- subroutine int_bounds1(total_ints,lower_bound,upper_bound) implicit none include 'DIMENSIONS' @@ -1566,6 +1593,47 @@ c--------------------------------------------------------------------------- endif return end +c--------------------------------------------------------------------------- + subroutine int_partition8(int_index,lower_index,upper_index,atom, + & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*) + implicit none + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + integer*8 int_index,lower_index,upper_index + integer atom,at_start,at_end, + & first_atom,last_atom,int_gr,jat_start,jat_end,int_index_old + logical lprn + lprn=.false. + if (lprn) write (iout,*) 'int_index=',int_index + int_index_old=int_index + int_index=int_index+last_atom-first_atom+1 + if (lprn) + & write (iout,*) 'int_index=',int_index, + & ' int_index_old',int_index_old, + & ' lower_index=',lower_index, + & ' upper_index=',upper_index, + & ' atom=',atom,' first_atom=',first_atom, + & ' last_atom=',last_atom + if (int_index.ge.lower_index) then + int_gr=int_gr+1 + if (at_start.eq.0) then + at_start=atom + jat_start=first_atom-1+lower_index-int_index_old + else + jat_start=first_atom + endif + if (lprn) write (iout,*) 'jat_start',jat_start + if (int_index.ge.upper_index) then + at_end=atom + jat_end=first_atom-1+upper_index-int_index_old + return1 + else + jat_end=last_atom + endif + if (lprn) write (iout,*) 'jat_end',jat_end + endif + return + end #endif c------------------------------------------------------------------------------ subroutine hpb_partition diff --git a/source/unres/src-HCD-5D/lagrangian_lesyng.F b/source/unres/src-HCD-5D/lagrangian_lesyng.F index 1180645..6dd113b 100644 --- a/source/unres/src-HCD-5D/lagrangian_lesyng.F +++ b/source/unres/src-HCD-5D/lagrangian_lesyng.F @@ -8,7 +8,7 @@ c------------------------------------------------------------------------- include 'DIMENSIONS' #ifdef MPI include 'mpif.h' - integer time00 + double precision time00 #endif include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -55,7 +55,7 @@ c------------------------------------------------------------------------- #endif #ifdef FIVEDIAG call grad_transform - d_a=0.0d0 + d_a(:,:2*nres)=0.0d0 if (lprn) then write (iout,*) "Potential forces backbone" do i=1,nres @@ -459,9 +459,9 @@ c write (iout,*) "i",i," itype",itype(i),ntyp1 endif enddo enddo - DMorig=DM - DU1orig=DU1 - DU2orig=DU2 + DMorig(:2*nres)=DM(:2*nres) + DU1orig(:2*nres)=DU1(:2*nres) + DU2orig(:2*nres)=DU2(:2*nres) if (gmatout) then write (iout,*)"The upper part of the five-diagonal inertia matrix" endif @@ -886,6 +886,9 @@ c--------------------------------------------------------------------------- c--------------------------------------------------------------------------- subroutine fivediaginv_mult(ndim,forces,d_a_vec) implicit none +#ifdef MPI + include 'mpif.h' +#endif include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' @@ -896,7 +899,12 @@ c--------------------------------------------------------------------------- double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim), & xsolv(ndim),d_a_vec(6*nres) integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev - accel=0.0d0 +#ifdef TIMING + include 'COMMON.TIME1' + double precision time01 + time01=MPI_Wtime() +#endif + accel(:,:2*nres)=0.0d0 do j=1,3 Compute accelerations in Calpha and SC do ichain=1,nchain @@ -996,6 +1004,9 @@ C Convert d_a to virtual-bon-vector basis ind=ind+3 endif enddo +#ifdef TIMING + time_ginvmult=time_ginvmult+MPI_Wtime()-time01 +#endif #ifdef DEBUG write (iout,*) "d_a_vec" write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside)) diff --git a/source/unres/src-HCD-5D/make_xx_list.F b/source/unres/src-HCD-5D/make_xx_list.F index 480aeb2..a1f6b45 100644 --- a/source/unres/src-HCD-5D/make_xx_list.F +++ b/source/unres/src-HCD-5D/make_xx_list.F @@ -22,6 +22,7 @@ ! print *,"START make_SC" #ifdef DEBUG write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res + write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e #endif r_buff_list=5.0d0 ilist_sc=0 @@ -100,9 +101,9 @@ #ifdef MPI #ifdef DEBUG write (iout,*) "before MPIREDUCE",ilist_sc - do i=1,ilist_sc - write (iout,*) i,contlisti(i),contlistj(i) - enddo +c do i=1,ilist_sc +c write (iout,*) i,contlisti(i),contlistj(i) +c enddo #endif if (nfgtasks.gt.1)then @@ -171,10 +172,226 @@ c write (iout,*) "SCSC after bcast ierr",ierr write (iout,*) i,newcontlisti(i),newcontlistj(i) enddo #endif - call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end) + call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end) +#ifdef DEBUG + write (iout,*) "g_listscsc_start",g_listscsc_start, + & "g_listscsc_end",g_listscsc_end return +#endif end -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine make_SCSC_inter_list_RESPA + implicit none + include "DIMENSIONS" +#ifdef MPI + include 'mpif.h' + include "COMMON.SETUP" +#endif + include "COMMON.CONTROL" + include "COMMON.CHAIN" + include "COMMON.INTERACT" + include "COMMON.SPLITELE" + include "COMMON.IOUNITS" + double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp + double precision dist_init, dist_temp,r_buff_list + integer contlisti_long(maxint_res*maxres), + & contlisti_short(maxint_res*maxres), + & contlistj_long(maxint_res*maxres), + & contlistj_short(maxint_res*maxres) +! integer :: newcontlisti(200*nres),newcontlistj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint, + & ilist_sc_long,g_ilist_sc_long,ilist_sc_short,g_ilist_sc_short + integer displ(0:max_fg_procs),i_ilist_sc_long(0:max_fg_procs), + & i_ilist_sc_short(0:max_fg_procs),ierr + logical lprn /.false./ + double precision boxshift + double precision d_scale,r_respa_buf +! print *,"START make_SC" +#ifdef DEBUG + write (iout,*) "make_SCSC_inter_list maxint_res",maxint_res + write (iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e +#endif + r_buff_list=5.0d0 + r_respa_buf=rlamb + ilist_sc_long=0 + ilist_sc_short=0 + do i=iatsc_s,iatsc_e + itypi=iabs(itype(i)) + if (itypi.eq.ntyp1) cycle + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + call to_box(xi,yi,zi) + do iint=1,nint_gr(i) + do j=istart(i,iint),iend(i,iint) + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + dist_init=dsqrt(xj*xj+yj*yj+zj*zj) +! r_buff_list is a read value for a buffer + if (dist_init.le.(r_cut_int+r_buff_list)) then +! Here the list is created + d_scale=dist_init/sigmaii(itypi,itypj) + if (d_scale.le.r_cut_respa+r_respa_buf) then + ilist_sc_short=ilist_sc_short+1 + contlisti_short(ilist_sc_short)=i + contlistj_short(ilist_sc_short)=j + endif + if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then + ilist_sc_long=ilist_sc_long+1 +! this can be substituted by cantor and anti-cantor + contlisti_long(ilist_sc_long)=i + contlistj_long(ilist_sc_long)=j + endif + endif + enddo + enddo + enddo +#ifdef MPI +#ifdef DEBUG + write (iout,*) "before MPIREDUCE ilist_sc_long",ilist_sc_long +c do i=1,ilist_sc_long +c write (iout,*) i,contlisti_long(i),contlistj_long(i) +c enddo + write (iout,*) "before MPIREDUCE ilist_sc_short",ilist_sc_short +c do i=1,ilist_sc_short +c write (iout,*) i,contlisti_short(i),contlistj_short(i) +c enddo +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_sc_long,g_ilist_sc_long,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) + call MPI_Reduce(ilist_sc_short,g_ilist_sc_short,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +c write (iout,*) "SCSC after reduce ierr",ierr + if (fg_rank.eq.0.and.(g_ilist_sc_long.gt.maxres*maxint_res .or. + & g_ilist_sc_short.gt.maxres*maxint_res)) then + if ((me.eq.king.or.out1file).and.energy_dec) then + write (iout,*) "Too many SCSC interactions", + & g_ilist_sc_long,g_ilist_sc_short, + & " only",maxres*maxint_res," allowed." + write (iout,*) "Specify a smaller r_cut_int and resubmit" + call flush(iout) + endif + write (*,*) "Processor:",me,": Too many SCSC interactions", + & g_ilist_sc_long+g_ilist_sc_short," only", + & maxres*maxint_res," allowed." + write (*,*) "Specify a smaller r_cut_int and resubmit" + call MPI_Abort(MPI_COMM_WORLD,ierr) + endif +c write(iout,*) "before bcast",g_ilist_sc_long + call MPI_Gather(ilist_sc_long,1,MPI_INTEGER, + & i_ilist_sc_long,1,MPI_INTEGER,king,FG_COMM,IERR) +c write (iout,*) "SCSC after gather ierr",ierr + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_sc_long(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlisti_long,ilist_sc_long,MPI_INTEGER, + & newcontlisti_long,i_ilist_sc_long,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCSC after gatherv ierr",ierr + call MPI_Gatherv(contlistj_long,ilist_sc_long,MPI_INTEGER, + & newcontlistj_long,i_ilist_sc_long,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_sc_long,1,MPI_INT,king,FG_COMM,IERR) +c write (iout,*) "SCSC bcast reduce ierr",ierr +! write(iout,*) "before bcast",g_ilist_sc_long +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + call MPI_Bcast(newcontlisti_long,g_ilist_sc_long,MPI_INT,king, + & FG_COMM,IERR) +c write (iout,*) "SCSC bcast reduce ierr",ierr + call MPI_Bcast(newcontlistj_long,g_ilist_sc_long,MPI_INT,king, + & FG_COMM,IERR) +c write (iout,*) "SCSC after bcast ierr",ierr +! write(iout,*) "before gather",displ(0),displ(1) +c write(iout,*) "before bcast",g_ilist_sc_short + call MPI_Gather(ilist_sc_short,1,MPI_INTEGER, + & i_ilist_sc_short,1,MPI_INTEGER,king,FG_COMM,IERR) +c write (iout,*) "SCSC after gather ierr",ierr + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_sc_short(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlisti_short,ilist_sc_short,MPI_INTEGER, + & newcontlisti_short,i_ilist_sc_short,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCSC after gatherv ierr",ierr + call MPI_Gatherv(contlistj_short,ilist_sc_short,MPI_INTEGER, + & newcontlistj_short,i_ilist_sc_short,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_sc_short,1,MPI_INT,king,FG_COMM,IERR) +c write (iout,*) "SCSC bcast reduce ierr",ierr +! write(iout,*) "before bcast",g_ilist_sc_short +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + call MPI_Bcast(newcontlisti_short,g_ilist_sc_short,MPI_INT,king, + & FG_COMM,IERR) +c write (iout,*) "SCSC bcast reduce ierr",ierr + call MPI_Bcast(newcontlistj_short,g_ilist_sc_short,MPI_INT,king, + ^ FG_COMM,IERR) +c write (iout,*) "SCSC after bcast ierr",ierr + else +#endif + g_ilist_sc_long=ilist_sc_long + + do i=1,ilist_sc_long + newcontlisti_long(i)=contlisti_long(i) + newcontlistj_long(i)=contlistj_long(i) + enddo + + g_ilist_sc_short=ilist_sc_short + + do i=1,ilist_sc_short + newcontlisti_short(i)=contlisti_short(i) + newcontlistj_short(i)=contlistj_short(i) + enddo +#ifdef MPI + endif +#endif + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + & write (iout,'(a30,2i10,a,2i4)') + & "Number of long- and short-range SC-SC interactions", + & g_ilist_sc_long,g_ilist_sc_short," per residue on average", + & g_ilist_sc_long/nres,g_ilist_sc_short/nres +#ifdef DEBUG + write (iout,*) + & "make_SCSC_inter_list: g_ilist_sc_long after GATHERV", + & g_ilist_sc_long + write (iout,*) "List of long-range SCSC interactions" + do i=1,g_ilist_sc_long + write (iout,*) i,newcontlisti_long(i),newcontlistj_long(i) + enddo + write (iout,*) + & "make_SCSC_inter_list: g_ilist_sc_short after GATHERV", + & g_ilist_sc_short + write (iout,*) "List of short-range SCSC interactions" + do i=1,g_ilist_sc_short + write (iout,*) i,newcontlisti_short(i),newcontlistj_short(i) + enddo +#endif + call int_bounds(g_ilist_sc_long,g_listscsc_start_long, + & g_listscsc_end_long) + call int_bounds(g_ilist_sc_short,g_listscsc_start_short, + & g_listscsc_end_short) +#ifdef DEBUG + write (iout,*) "g_list_sc_start",g_listscsc_start_long, + & "g_list_sc_end",g_listscsc_end_long + write (iout,*)"g_list_sc_start_short",g_listscsc_start_short, + & "g_list_sc_end_short",g_listscsc_end_short +#endif + return + end +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine make_SCp_inter_list implicit none include "DIMENSIONS" @@ -204,8 +421,8 @@ c & contlistscpj_f(2*maxint_res*maxres) write (iout,*) "make_SCp_inter_list maxint_res",maxint_res #endif r_buff_list=5.0 - ilist_scp=0 - ilist_scp_first=0 + ilist_scp=0 + ilist_scp_first=0 do i=iatscp_s,iatscp_e if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle xi=0.5D0*(c(1,i)+c(1,i+1)) @@ -377,11 +594,234 @@ c write (iout,*) "SCp bcast reduce ierr",ierr ! enddo #endif call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end) +#ifdef DEBUG + write (iout,*) "g_listscp_start",g_listscp_start, + & "g_listscp_end",g_listscp_end +#endif + return + end +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine make_SCp_inter_list_RESPA + implicit none + include "DIMENSIONS" +#ifdef MPI + include 'mpif.h' + include "COMMON.SETUP" +#endif + include "COMMON.CONTROL" + include "COMMON.CHAIN" + include "COMMON.INTERACT" + include "COMMON.SPLITELE" + include "COMMON.IOUNITS" + double precision xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp + double precision dist_init, dist_temp,r_buff_list + integer contlistscpi_long(2*maxint_res*maxres), + & contlistscpi_short(2*maxint_res*maxres), + & contlistscpj_long(2*maxint_res*maxres), + & contlistscpj_short(2*maxint_res*maxres) +! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres) + integer i,j,iteli,itypj,subchap,xshift,yshift,zshift,iint, + & ilist_scp_long,ilist_scp_short,g_ilist_scp_long,g_ilist_scp_short + integer displ(0:max_fg_procs),i_ilist_scp_long(0:max_fg_procs), + & i_ilist_scp_short(0:max_fg_procs),ierr +c integer contlistscpi_f(2*maxint_res*maxres), +c & contlistscpj_f(2*maxint_res*maxres) + double precision boxshift + double precision d_scale,r_respa_buf +! print *,"START make_SC" +#ifdef DEBUG + write (iout,*) "make_SCp_inter_list maxint_res",maxint_res +#endif + r_buff_list=5.0 + r_respa_buf=rlamb + ilist_scp_long=0 + ilist_scp_short=0 + do i=iatscp_s,iatscp_e + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + xi=0.5D0*(c(1,i)+c(1,i+1)) + yi=0.5D0*(c(2,i)+c(2,i+1)) + zi=0.5D0*(c(3,i)+c(3,i+1)) + call to_box(xi,yi,zi) + iteli=itel(i) + do iint=1,nscp_gr(i) + do j=iscpstart(i,iint),iscpend(i,iint) + itypj=iabs(itype(j)) + if (itypj.eq.ntyp1) cycle +! Uncomment following three lines for SC-p interactions +! xj=c(1,nres+j)-xi +! yj=c(2,nres+j)-yi +! zj=c(3,nres+j)-zi +! Uncomment following three lines for Ca-p interactions +! xj=c(1,j)-xi +! yj=c(2,j)-yi +! zj=c(3,j)-zi + xj=c(1,j) + yj=c(2,j) + zj=c(3,j) + call to_box(xj,yj,zj) + xj=boxshift(xj-xi,boxxsize) + yj=boxshift(yj-yi,boxysize) + zj=boxshift(zj-zi,boxzsize) + dist_init=dsqrt(xj*xj+yj*yj+zj*zj) +! r_buff_list is a read value for a buffer + if (dist_init.le.(r_cut_int+r_buff_list)) then + d_scale=dist_init/rscp(itypj,iteli) + if (d_scale.le.r_cut_respa+r_respa_buf) then +! Here the list is created + ilist_scp_short=ilist_scp_short+1 + contlistscpi_short(ilist_scp_short)=i + contlistscpj_short(ilist_scp_short)=j + endif + if (d_scale.gt.r_cut_respa-rlamb-r_respa_buf) then +! this can be substituted by cantor and anti-cantor + ilist_scp_long=ilist_scp_long+1 + contlistscpi_long(ilist_scp_long)=i + contlistscpj_long(ilist_scp_long)=j + endif + endif + enddo + enddo + enddo +#ifdef MPI +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_scp_long,ilist_scp_short + write (iout,*) "Long-range scp interaction list" + do i=1,ilist_scp_long + write (iout,*) i,contlistscpi_long(i),contlistscpj_long(i) + enddo + write (iout,*) "Short-range scp interaction list" + do i=1,ilist_scp_short + write (iout,*) i,contlistscpi_short(i),contlistscpj_short(i) + enddo +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_scp_long,g_ilist_scp_long,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) + call MPI_Reduce(ilist_scp_short,g_ilist_scp_short,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +c write (iout,*) "SCp after reduce ierr",ierr + if (fg_rank.eq.0.and.(g_ilist_scp_long.gt. + & 2*maxres*maxint_res .or. g_ilist_scp_short.gt. + & 2*maxres*maxint_res)) then + if ((me.eq.king.or.out1file).and.energy_dec) then + write (iout,*) "Too many SCp interactions", + & g_ilist_scp_long+g_ilist_scp_short," only", + & 2*maxres*maxint_res," allowed." + write (iout,*) "Specify a smaller r_cut_int and resubmit" + call flush(iout) + endif + write (*,*) "Processor:",me,": Too many SCp interactions", + & g_ilist_scp_long+g_ilist_scp_short," only", + & 2*maxres*maxint_res," allowed." + write (*,*) "Specify a smaller r_cut_int and resubmit" + call MPI_Abort(MPI_COMM_WORLD,ierr) + endif +c write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_scp_long,1,MPI_INTEGER, + & i_ilist_scp_long,1,MPI_INTEGER,king,FG_COMM,IERR) + call MPI_Gather(ilist_scp_short,1,MPI_INTEGER, + & i_ilist_scp_short,1,MPI_INTEGER,king,FG_COMM,IERR) +c write (iout,*) "SCp after gather ierr",ierr + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_scp_long(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistscpi_long,ilist_scp_long,MPI_INTEGER, + & newcontlistscpi_long,i_ilist_scp_long,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after gatherv ierr",ierr + call MPI_Gatherv(contlistscpj_long,ilist_scp_long,MPI_INTEGER, + & newcontlistscpj_long,i_ilist_scp_long,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after gatherv ierr",ierr + call MPI_Bcast(g_ilist_scp_long,1,MPI_INT,king,FG_COMM,IERR) +c write (iout,*) "SCp after bcast ierr",ierr +! write(iout,*) "before bcast",g_ilist_sc +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + call MPI_Bcast(newcontlistscpi_long,g_ilist_scp_long,MPI_INT, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after bcast ierr",ierr + call MPI_Bcast(newcontlistscpj_long,g_ilist_scp_long,MPI_INT, + & king,FG_COMM,IERR) +c write (iout,*) "SCp bcast reduce ierr",ierr +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_scp_short(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistscpi_short,ilist_scp_short,MPI_INTEGER, + & newcontlistscpi_short,i_ilist_scp_short,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after gatherv ierr",ierr + call MPI_Gatherv(contlistscpj_short,ilist_scp_short,MPI_INTEGER, + & newcontlistscpj_short,i_ilist_scp_short,displ,MPI_INTEGER, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after gatherv ierr",ierr + call MPI_Bcast(g_ilist_scp_short,1,MPI_INT,king,FG_COMM,IERR) + call MPI_Bcast(newcontlistscpi_short,g_ilist_scp_short,MPI_INT, + & king,FG_COMM,IERR) +c write (iout,*) "SCp after bcast ierr",ierr + call MPI_Bcast(newcontlistscpj_short,g_ilist_scp_short,MPI_INT, + & king,FG_COMM,IERR) + else +#endif + g_ilist_scp_long=ilist_scp_long + + do i=1,ilist_scp_long + newcontlistscpi_long(i)=contlistscpi_long(i) + newcontlistscpj_long(i)=contlistscpj_long(i) + enddo + g_ilist_scp_short=ilist_scp_short + + do i=1,ilist_scp_short + newcontlistscpi_short(i)=contlistscpi_short(i) + newcontlistscpj_short(i)=contlistscpj_short(i) + enddo +#ifdef MPI + endif +#endif + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + &then + write (iout,'(a30,i10,a,i4)') + & "Number of long-range SC-p interactions", + & g_ilist_scp_long," per residue on average",g_ilist_scp_long/nres + write (iout,'(a30,i10,a,i4)') + & "Number of short-range SC-p interactions", + &g_ilist_scp_short," per residue on average",g_ilist_scp_short/nres + endif +#ifdef DEBUG + write (iout,*) "make_SCp_inter_list: after GATHERV long-range", + & g_ilist_scp_long + do i=1,g_ilist_scp_long + write (iout,*) i,newcontlistscpi_long(i),newcontlistscpj_long(i) + enddo + write (iout,*) "make_SCp_inter_list: after GATHERV short-range", + & g_ilist_scp_short + do i=1,g_ilist_scp_short + write (iout,*) i,newcontlistscpi_short(i), + & newcontlistscpj_short(i) + enddo +#endif + call int_bounds(g_ilist_scp_long,g_listscp_start_long, + & g_listscp_end_long) + call int_bounds(g_ilist_scp_short,g_listscp_start_short, + & g_listscp_end_short) + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + &then + write (iout,*) "g_listscp_start",g_listscp_start_long, + & "g_listscp_end",g_listscp_end_long + write (iout,*)"g_listscp_start_short",g_listscp_start_short, + & "g_listscp_end_short",g_listscp_end_short + endif return end !----------------------------------------------------------------------------- - subroutine make_pp_vdw_inter_list + subroutine make_pp_vdw_inter_list_RESPA implicit none include "DIMENSIONS" #ifdef MPI @@ -398,164 +838,148 @@ c write (iout,*) "SCp bcast reduce ierr",ierr double precision xmedj,ymedj,zmedj double precision dist_init, dist_temp,r_buff_list,dxi,dyi,dzi, & xmedi,ymedi,zmedi - double precision dx_normi,dy_normi,dz_normi,dxj,dyj,dzj, - & dx_normj,dy_normj,dz_normj - integer contlistpp_vdwi(maxint_res*maxres), - & contlistpp_vdwj(maxint_res*maxres) -! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres) + double precision dxj,dyj,dzj + integer contlistpp_vdwi_short(maxint_res*maxres), + & contlistpp_vdwj_short(maxint_res*maxres) integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint, - & ilist_pp_vdw,g_ilist_pp_vdw - integer displ(0:max_fg_procs),i_ilist_pp_vdw(0:max_fg_procs),ierr + & ilist_pp_vdw_short,g_ilist_pp_vdw_short + integer displ(0:max_fg_procs), + & i_ilist_pp_vdw_short(0:max_fg_procs),ierr ! print *,"START make_SC" + double precision boxshift + double precision d_scale,r_respa_buf #ifdef DEBUG write (iout,*) "make_pp_vdw_inter_list" #endif - ilist_pp_vdw=0 + ilist_pp_vdw_short=0 r_buff_list=5.0 + r_respa_buf=rlamb do i=iatel_s_vdw,iatel_e_vdw if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) - dx_normi=dc_norm(1,i) - dy_normi=dc_norm(2,i) - dz_normi=dc_norm(3,i) xmedi=c(1,i)+0.5d0*dxi ymedi=c(2,i)+0.5d0*dyi zmedi=c(3,i)+0.5d0*dzi - xmedi=dmod(xmedi,boxxsize) - if (xmedi.lt.0) xmedi=xmedi+boxxsize - ymedi=dmod(ymedi,boxysize) - if (ymedi.lt.0) ymedi=ymedi+boxysize - zmedi=dmod(zmedi,boxzsize) - if (zmedi.lt.0) zmedi=zmedi+boxzsize + call to_box(xmedi,ymedi,zmedi) do j=ielstart_vdw(i),ielend_vdw(i) ! write (iout,*) i,j,itype(i),itype(j) if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle - ! 1,j) dxj=dc(1,j) dyj=dc(2,j) dzj=dc(3,j) - dx_normj=dc_norm(1,j) - dy_normj=dc_norm(2,j) - dz_normj=dc_norm(3,j) ! xj=c(1,j)+0.5D0*dxj-xmedi ! yj=c(2,j)+0.5D0*dyj-ymedi ! zj=c(3,j)+0.5D0*dzj-zmedi xj=c(1,j)+0.5D0*dxj yj=c(2,j)+0.5D0*dyj zj=c(3,j)+0.5D0*dzj - xj=mod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=mod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=mod(zj,boxzsize) - if (zj.lt.0) zj=zj+boxzsize - - dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - endif - enddo - enddo - enddo + call to_box(xj,yj,zj) + xj=boxshift(xj-xmedi,boxxsize) + yj=boxshift(yj-ymedi,boxysize) + zj=boxshift(zj-zmedi,boxzsize) + dist_init=dsqrt(xj*xj+yj*yj+zj*zj) - if (dsqrt(dist_init).le.(r_cut_int+r_buff_list)) then + if (dist_init.le.(r_cut_int+r_buff_list)) then + d_scale=dist_init/rpp(itel(i),itel(j)) + if (d_scale.le.r_cut_respa+r_respa_buf) then ! Here the list is created - ilist_pp_vdw=ilist_pp_vdw+1 + ilist_pp_vdw_short=ilist_pp_vdw_short+1 ! this can be substituted by cantor and anti-cantor - contlistpp_vdwi(ilist_pp_vdw)=i - contlistpp_vdwj(ilist_pp_vdw)=j + contlistpp_vdwi_short(ilist_pp_vdw_short)=i + contlistpp_vdwj_short(ilist_pp_vdw_short)=j + endif endif - enddo - enddo + enddo + enddo ! enddo #ifdef MPI #ifdef DEBUG - write (iout,*) "before MPIREDUCE",ilist_pp_vdw - do i=1,ilist_pp_vdw - write (iout,*) i,contlistpp_vdwi(i),contlistpp_vdwj(i) + write (iout,*) "before MPIREDUCE longrange",ilist_pp_vdw_long + do i=1,ilist_pp_vdw_long + write (iout,*) i,contlistpp_vdwi_long(i),contlistpp_vdwj_long(i) + enddo + write (iout,*) "before MPIREDUCE shortrange",ilist_pp_vdw_short + do i=1,ilist_pp_vdw_short + write (iout,*) i,contlistpp_vdwi_short(i), + & contlistpp_vdwj_short(i) enddo #endif if (nfgtasks.gt.1)then - call MPI_Reduce(ilist_pp_vdw,g_ilist_pp_vdw,1, + call MPI_Reduce(ilist_pp_vdw_short,g_ilist_pp_vdw_short,1, & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) - if (fg_rank.eq.0.and.g_ilist_pp_vdw.gt.maxres*maxint_res) then + if (fg_rank.eq.0.and.g_ilist_pp_vdw_short.gt.maxres*maxint_res) + & then if ((me.eq.king.or.out1file).and.energy_dec) then write (iout,*) "Too many pp VDW interactions", - & g_ilist_pp_vdw," only",maxres*maxint_res," allowed." + & g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed." write (iout,*) "Specify a smaller r_cut_int and resubmit" call flush(iout) endif write (*,*) "Processor:",me,": Too many pp VDW interactions", - & g_ilist_pp_vdw," only",maxres*maxint_res," allowed." + & g_ilist_pp_vdw_short," only",maxres*maxint_res," allowed." write (8,*) "Specify a smaller r_cut_int and resubmit" call MPI_Abort(MPI_COMM_WORLD,ierr) endif ! write(iout,*) "before bcast",g_ilist_sc - call MPI_Gather(ilist_pp_vdw,1,MPI_INTEGER, - & i_ilist_pp_vdw,1,MPI_INTEGER,king,FG_COMM,IERR) + call MPI_Gather(ilist_pp_vdw_short,1,MPI_INTEGER, + & i_ilist_pp_vdw_short,1,MPI_INTEGER,king,FG_COMM,IERR) displ(0)=0 do i=1,nfgtasks-1,1 - displ(i)=i_ilist_pp_vdw(i-1)+displ(i-1) + displ(i)=i_ilist_pp_vdw_short(i-1)+displ(i-1) enddo ! write(iout,*) "before gather",displ(0),displ(1) - call MPI_Gatherv(contlistpp_vdwi,ilist_pp_vdw,MPI_INTEGER, - & newcontlistpp_vdwi,i_ilist_pp_vdw,displ,MPI_INTEGER, - & king,FG_COMM,IERR) - call MPI_Gatherv(contlistpp_vdwj,ilist_pp_vdw,MPI_INTEGER, - & newcontlistpp_vdwj,i_ilist_pp_vdw,displ,MPI_INTEGER, - & king,FG_COMM,IERR) - call MPI_Bcast(g_ilist_pp_vdw,1,MPI_INT,king,FG_COMM,IERR) + call MPI_Gatherv(contlistpp_vdwi_short,ilist_pp_vdw_short, + & MPI_INTEGER,newcontlistpp_vdwi_short,i_ilist_pp_vdw_short,displ, + & MPI_INTEGER,king,FG_COMM,IERR) + call MPI_Gatherv(contlistpp_vdwj_short,ilist_pp_vdw_short, + & MPI_INTEGER,newcontlistpp_vdwj_short,i_ilist_pp_vdw_short,displ, + & MPI_INTEGER,king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_pp_vdw_short,1,MPI_INT,king,FG_COMM,IERR) ! write(iout,*) "before bcast",g_ilist_sc ! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) - call MPI_Bcast(newcontlistpp_vdwi,g_ilist_pp_vdw,MPI_INT,king, - & FG_COMM,IERR) - call MPI_Bcast(newcontlistpp_vdwj,g_ilist_pp_vdw,MPI_INT,king, - & FG_COMM,IERR) - + call MPI_Bcast(newcontlistpp_vdwi_short,g_ilist_pp_vdw_short, + & MPI_INT,king,FG_COMM,IERR) + call MPI_Bcast(newcontlistpp_vdwj_short,g_ilist_pp_vdw_short, + & MPI_INT,king,FG_COMM,IERR) ! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) - - else + else #endif - g_ilist_pp_vdw=ilist_pp_vdw + g_ilist_pp_vdw_short=ilist_pp_vdw_short - do i=1,ilist_pp_vdw - newcontlistpp_vdwi(i)=contlistpp_vdwi(i) - newcontlistpp_vdwj(i)=contlistpp_vdwj(i) + do i=1,ilist_pp_vdw_short + newcontlistpp_vdwi_short(i)=contlistpp_vdwi_short(i) + newcontlistpp_vdwj_short(i)=contlistpp_vdwj_short(i) enddo #ifdef MPI - endif + endif #endif - call int_bounds(g_ilist_pp_vdw,g_listpp_vdw_start, - & g_listpp_vdw_end) - if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) - &write (iout,'(a30,i10,a,i4)') "Number of p-p VDW interactions", - & g_ilist_pp_vdw," per residue on average",g_ilist_pp_vdw/nres + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + &then + write (iout,*) "Number of short-range p-p VDW interactions", + & g_ilist_pp_vdw_short," per residue on average", + & g_ilist_pp_vdw_short/nres + endif #ifdef DEBUG - write (iout,*) "g_listpp_vdw_start",g_listpp_vdw_start, - & "g_listpp_vdw_end",g_listpp_vdw_end + write (iout,*) "Short-range pp_vdw" write (iout,*) "make_pp_vdw_inter_list: after GATHERV", - & g_ilist_pp_vdw - do i=1,g_ilist_pp_vdw - write (iout,*) i,newcontlistpp_vdwi(i),newcontlistpp_vdwj(i) + & g_ilist_pp_vdw_short + do i=1,g_ilist_pp_vdw_short + write (iout,*) i,newcontlistpp_vdwi_short(i), + & newcontlistpp_vdwj_short(i) enddo #endif + call int_bounds(g_ilist_pp_vdw_short,g_listpp_vdw_start_short, + & g_listpp_vdw_end_short) + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + &then + write (iout,*)"g_listpp_vdw_start_short", + & g_listpp_vdw_start_short, + & "g_listpp_vdw_end_short",g_listpp_vdw_end_short + endif return end !----------------------------------------------------------------------------- @@ -716,18 +1140,17 @@ c write (iout,*) "After bcast ierr",ierr ! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) - else + else #endif g_ilist_pp=ilist_pp do i=1,ilist_pp - newcontlistppi(i)=contlistppi(i) - newcontlistppj(i)=contlistppj(i) + newcontlistppi(i)=contlistppi(i) + newcontlistppj(i)=contlistppj(i) enddo #ifdef MPI - endif + endif #endif - call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end) if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) & write (iout,'(a30,i10,a,i4)') "Number of p-p interactions", & g_ilist_pp," per residue on average",g_ilist_pp/nres @@ -737,5 +1160,11 @@ c write (iout,*) "After bcast ierr",ierr write (iout,*) i,newcontlistppi(i),newcontlistppj(i) enddo #endif + call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end) + if (fg_rank.eq.0 .and. (me.eq.king.or.out1file) .and. energy_dec) + &then + write (iout,*) "g_listpp_start",g_listpp_start, + & "g_listpp_end",g_listpp_end + endif return end diff --git a/source/unres/src-HCD-5D/minimize_p.F b/source/unres/src-HCD-5D/minimize_p.F index 41a1a27..0b16d55 100644 --- a/source/unres/src-HCD-5D/minimize_p.F +++ b/source/unres/src-HCD-5D/minimize_p.F @@ -197,11 +197,11 @@ c---------------------------------------------------------------------------- #endif include 'COMMON.TIME1' double precision z(maxres6),d_a_tmp(maxres6) - double precision edum(0:n_ene),time_order(0:10) + double precision edum(0:n_ene),time_order(0:11) c double precision Gcopy(maxres2,maxres2) c common /przechowalnia/ Gcopy integer icall /0/ - integer i,j,iorder + integer i,j,iorder,ioverlap(maxres),ioverlap_last C Workers wait for variables and NF, and NFL from the boss iorder=0 do while (iorder.ge.0) @@ -301,6 +301,8 @@ c call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block) #endif else if (iorder.eq.10) then call setup_fricmat + else if (iorder.eq.11) then + call overlap_sc_list(ioverlap,ioverlap_last,.false.) endif enddo write (*,*) 'Processor',fg_rank,' CG group',kolor, diff --git a/source/unres/src-HCD-5D/orig_frame_chain.F b/source/unres/src-HCD-5D/orig_frame_chain.F new file mode 100644 index 0000000..07053ce --- /dev/null +++ b/source/unres/src-HCD-5D/orig_frame_chain.F @@ -0,0 +1,85 @@ + subroutine orig_frame_chain(istart) +C +C Define the origin and orientation of the coordinate system starting +C at residue istart and locate sites istart+1 and istart+2. The +C coordinates of site istart and the respective dc and dc_norm must be +C pre-defined +C + implicit none + integer i,j,istart,ichain + double precision cost,sint,cosg,sing,aux + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.LOCAL' + include 'COMMON.GEO' + include 'COMMON.VAR' + cost=dc_norm(1,istart) + aux=dsqrt(dc_norm(2,istart)**2+dc_norm(3,istart)**2) + cosg=dc_norm(2,istart)/aux + sing=dc_norm(3,istart)/aux + cost=dcos(theta(istart+2)) + sint=dsin(theta(istart+2)) + prod(1,1,istart)=cost + prod(1,2,istart)=-sint + prod(1,3,istart)=0.0d0 + prod(2,1,istart)=sint*cosg + prod(2,2,istart)=cost*cosg + prod(2,3,istart)=-sing + prod(3,1,istart)=sint*sing + prod(3,2,istart)=cost*sing + prod(3,3,istart)=cosg + aux=prod(1,1,istart)*(prod(2,2,istart)*prod(3,3,istart) + & -prod(3,2,istart)*prod(2,3,istart)) + & -prod(1,2,istart)*(prod(2,1,istart)*prod(3,3,istart) + & -prod(3,1,istart)*prod(2,3,istart)) + & +prod(1,3,istart)*(prod(2,1,istart)*prod(3,2,istart) + & -prod(3,1,istart)*prod(2,2,istart)) +c write (iout,*) "orig_frame_chain prod",istart +c do i=1,3 +c write(iout,'(i5,3f10.5)') i,(prod(i,j,istart),j=1,3) +c enddo +c write (iout,*) "orig_frame_chain: prod",istart," determinant",aux + t(1,1,istart)=cost + t(1,2,istart)=-sint + t(1,3,istart)= 0.0D0 + t(2,1,istart)=sint + t(2,2,istart)=cost + t(2,3,istart)= 0.0D0 + t(3,1,istart)= 0.0D0 + t(3,2,istart)= 0.0D0 + t(3,3,istart)= 1.0D0 + r(1,1,istart)= 1.0D0 + r(1,2,istart)= 0.0D0 + r(1,3,istart)= 0.0D0 + r(2,1,istart)= 0.0D0 + r(2,2,istart)= 1.0D0 + r(2,3,istart)= 0.0D0 + r(3,1,istart)= 0.0D0 + r(3,2,istart)= 0.0D0 + r(3,3,istart)= 1.0D0 + do i=1,3 + do j=1,3 + rt(i,j,istart)=t(i,j,istart) + enddo + enddo + call matmult(prod(1,1,istart),rt(1,1,istart),prod(1,1,istart+1)) +c aux=prod(1,1,istart+1)*(prod(2,2,istart+1)*prod(3,3,istart+1) +c & -prod(3,2,istart+1)*prod(2,3,istart+1)) +c & -prod(1,2,istart+1)*(prod(2,1,istart+1)*prod(3,3,istart+1) +c & -prod(3,1,istart+1)*prod(2,3,istart+1)) +c & +prod(1,3,istart+1)*(prod(2,1,istart+1)*prod(3,2,istart+1) +c & -prod(3,1,istart+1)*prod(2,2,istart+1)) +c write (iout,*) "orig_frame_chain prod",istart+1 +c do i=1,3 +c write(iout,'(i5,3f10.5)') i,(prod(i,j,istart+1),j=1,3) +c enddo +c write (iout,*)"orig_frame_chain: prod",istart+1," determinant",aux + do j=1,3 + dc_norm(j,istart+1)=prod(j,1,istart+1) + dc(j,istart+1)=vbld(istart+2)*prod(j,1,istart+1) + c(j,istart+2)=c(j,istart+1)+dc(j,istart+1) + enddo + call locate_side_chain(istart+1) + return + end diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F index 41fe7f6..12d011e 100644 --- a/source/unres/src-HCD-5D/readpdb-mult.F +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -64,9 +64,9 @@ c call flush(iout) iterter(ires_old-1)=1 itype(ires_old)=ntyp1 iterter(ires_old)=1 - ishift1=ishift1+1 +c ishift1=ishift1+1 ibeg=2 - write (iout,*) "Chain ended",ires,ishift,ires_old + write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) @@ -95,8 +95,8 @@ c write (iout,*) "! ",atom," !",ires read (card(18:20),'(a3)') res c write (iout,*) "ires",ires,ires-ishift+ishift1, c & " ires_old",ires_old -c write (iout,*) "ishift",ishift," ishift1",ishift1 -c write (iout,*) "IRES",ires-ishift+ishift1,ires_old +c write (iout,*) "ishift",ishift," ishift1",ishift1 +c write (iout,*) "IRES",ires-ishift+ishift1,ires_old if (ires-ishift+ishift1.ne.ires_old) then ! Calculate the CM of the preceding residue. ! if (ibeg.eq.0) call sccenter(ires,iii,sccor) @@ -115,6 +115,7 @@ c write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3) sccalc=.true. endif ! Start new residue. +c write (iout,*) "ibeg",ibeg if (res.eq.'Cl-' .or. res.eq.'Na+') then ires=ires_old cycle @@ -133,10 +134,13 @@ c write (iout,*) "BEG ires",ires else if (ibeg.eq.2) then ! Start a new chain ishift=-ires_old+ires-1 !!!!! - ishift1=ishift1-1 !!!!! -c write (iout,*) "New chain started",ires,ishift,ishift1,"!" +c ishift1=ishift1-1 !!!!! +c write (iout,*) "New chain started",ires,ires_old,ishift, +c & ishift1 ires=ires-ishift+ishift1 + write (iout,*) "New chain started ires",ires ires_old=ires +c ires=ires_old+1 ibeg=0 else ishift=ishift-(ires-ishift+ishift1-ires_old-1) @@ -159,7 +163,9 @@ c write (2,*) "ires",ires," res ",res!," ity"!,ity if (atom.eq.'CA' .or. atom.eq.'CH3' .or. & res.eq.'NHE'.and.atom(:2).eq.'HN') then read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) -! write (iout,*) "backbone ",atom + read(card(61:66),*) bfac(ires) +c write (iout,*) "backbone ",atom +c write (iout,*) ires,res,(c(j,ires),j=1,3) #ifdef DEBUG write (iout,'(i6,i3,2x,a,3f8.3)') & ires,itype(ires),res,(c(j,ires),j=1,3) @@ -181,7 +187,7 @@ c write (2,*) "iii",iii endif enddo 10 if(me.eq.king.or..not.out1file) - & write (iout,'(a,i5)') ' Nres: ',ires + & write (iout,'(a,i7)') ' Nres: ',ires c write (iout,*) "iii",iii C Calculate dummy residue coordinates inside the "chain" of a multichain C system @@ -229,7 +235,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else !unres_pdb do j=1,3 @@ -316,9 +322,9 @@ C Calculate internal coordinates. write (iout,'(/a)') & "Cartesian coordinates of the reference structure" write (iout,'(a,3(3x,a5),5x,3(3x,a5))') - & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" + & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)" do ires=1,nres - write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') + write (iout,'(a3,1x,i6,3f8.3,5x,3f8.3)') & restyp(itype(ires)),ires,(c(j,ires),j=1,3), & (c(j,ires+nres),j=1,3) enddo @@ -412,8 +418,9 @@ C and convert the peptide geometry into virtual-chain geometry. character*3 seq,res character*5 atom character*80 card - double precision sccor(3,20) + double precision sccor(3,50) integer rescode,iterter(maxres) + logical zero do i=1,maxres iterter(i)=0 enddo @@ -541,7 +548,7 @@ c write (iout,*) "sidechain ",atom endif enddo 10 if(me.eq.king.or..not.out1file) - & write (iout,'(a,i5)') ' Nres: ',ires + & write (iout,'(a,i7)') ' Nres: ',ires C Calculate dummy residue coordinates inside the "chain" of a multichain C system nres=ires @@ -581,7 +588,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,i)=c(j,i+1)-1.9d0*e2(j) + c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else !unres_pdb do j=1,3 @@ -616,7 +623,7 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,nres)=c(j,nres-1)-1.9d0*e2(j) + c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo else do j=1,3 @@ -648,7 +655,7 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue e2(3)=0.0d0 endif do j=1,3 - c(j,1)=c(j,2)-1.9d0*e2(j) + c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0) enddo else do j=1,3 @@ -676,6 +683,18 @@ C Calculate internal coordinates. & (c(j,ires+nres),j=1,3) enddo endif + zero=.false. + do ires=1,nres + zero=zero.or.itype(ires).eq.0 + enddo + if (zero) then + write (iout,'(2a)') "Gaps in PDB coordinates detected;", + & " look for ZERO in the control output above." + write (iout,'(2a)') "Repair the PDB file using MODELLER", + & " or other softwared and resubmit." + call flush(iout) + stop + endif C Calculate internal coordinates. call int_from_cart(.true.,out_template_coord) call sc_loc_geom(.false.) @@ -683,6 +702,7 @@ C Calculate internal coordinates. thetaref(i)=theta(i) phiref(i)=phi(i) enddo + dc(:,0)=c(:,1) do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 4fbc0f1..58b86cc 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -149,6 +149,7 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0 call readi(controlcard,'SYM',symetr,1) + call readi(controlcard,'PERMUT',npermut,1) call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes @@ -167,6 +168,8 @@ c write (iout,'(a,f10.1)') 'Time limit (min):',timlim c endif call readi(controlcard,'NZ_START',nz_start,0) call readi(controlcard,'NZ_END',nz_end,0) + call readi(controlcard,'NRAN_START',nran_start,0) + write (iout,*) "nran_start",nran_start c call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety @@ -481,6 +484,7 @@ C call readi(controlcard,"NSTEP",n_timestep,1000000) call readi(controlcard,"NTWE",ntwe,100) call readi(controlcard,"NTWX",ntwx,1000) + call readi(controlcard,"REST_FREQ",irest_freq,1000) call reada(controlcard,"DT",d_time,1.0d-1) call reada(controlcard,"DVMAX",dvmax,2.0d1) call reada(controlcard,"DAMAX",damax,1.0d1) @@ -556,6 +560,7 @@ c if performing umbrella sampling, fragments constrained are read from the frag & " A" write(iout,'(a60,i5)')"Frequency of updating interaction list", & imatupdate + write(iout,'(a60,i5)')"Restart writing frequency",irest_freq if (RESPA) then write (iout,'(2a,i4,a)') & "A-MTS algorithm used; initial time step for fast-varying", @@ -741,7 +746,7 @@ C integer ilen external ilen integer iperm,tperm - integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp,itemp + integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp double precision sumv C C Read PDB structure if applicable @@ -857,16 +862,19 @@ cd print *,'NNT=',NNT,' NCT=',NCT & chain_border(2,i),chain_border1(1,i),chain_border1(2,i) enddo call chain_symmetry(nchain,nres,itype,chain_border, - & chain_length,npermchain,tabpermchain) + & chain_length,npermchain,tabpermchain,nchain_group,nequiv, + & iequiv,mapchain) c do i=1,nres c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain), c & ii=1,npermchain) c enddo +#ifdef DEBUG write(iout,*) "residue permutations" do i=1,nres write(iout,*) i,(iperm(i,ii),ii=1,npermchain) enddo call flush(iout) +#endif if (itype(1).eq.ntyp1) nnt=2 if (itype(nres).eq.ntyp1) nct=nct-1 write (iout,*) "nnt",nnt," nct",nct @@ -971,25 +979,34 @@ C & write (iout,*) 'FTORS',ftors call reada(weightcard,"WDIHC",wdihc,0.591D0) write (iout,*) "Weight of dihedral angle restraints",wdihc read(inp,'(9x,3f7.3)') +c & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) write (iout,*) "The secprob array" do i=nnt,nct write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3) enddo ndih_constr=0 + iconstr_dih=0 do i=nnt+3,nct if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then ndih_constr=ndih_constr+1 idih_constr(ndih_constr)=i + iconstr_dih(i)=ndih_constr sumv=0.0d0 do j=1,3 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2) sumv=sumv+vpsipred(j,ndih_constr) enddo - do j=1,3 - vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv - enddo + if (sumv.gt.0.0d0) then + do j=1,3 + vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv + enddo + else + vpsipred(1,ndih_constr)=1.0d0 + vpsipred(2,ndih_constr)=0.0d0 + vpsipred(3,ndih_constr)=0.0d0 + endif phibound(1,ndih_constr)=phihel*deg2rad phibound(2,ndih_constr)=phibet*deg2rad sdihed(1,ndih_constr)=sigmahel*deg2rad @@ -1125,8 +1142,11 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb endif write (iout,*) "calling read_saxs_consrtr",nsaxs if (nsaxs.gt.0) call read_saxs_constr - +c write (iout,*) "After read_saxs_constr" +c call flush(iout) if (constr_homology.gt.0) then +c write (iout,*) "Calling read_constr_homology" +c call flush(iout) call read_constr_homology if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres @@ -1184,7 +1204,8 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb 335 continue unres_pdb=.false. nres_temp=nres - call readpdb +c call readpdb + call readpdb_template(nmodel_start+1) close(ipdbin) if (nres.ge.nres_temp) then nmodel_start=nmodel_start+1 @@ -1195,22 +1216,11 @@ c write (iout,*) "After read_dist_constr nhpb",nhpb enddo enddo else -c itemp=nres -c nres=nres_temp -c call gen_rand_conf(itemp,*115) -c nmodel_start=nmodel_start+1 -c do i=1,2*nres -c do j=1,3 -c chomo(j,i,nmodel_start)=c(j,i) -c enddo -c enddo -c goto 116 - 115 if (me.eq.king .or. .not. out1file) - & write (iout,'(a,2i5,1x,a)') + if (me.eq.king .or. .not. out1file) + & write (iout,'(a,2i7,1x,a)') & "Different number of residues",nres_temp,nres, & " model skipped." endif - 116 continue nres=nres_temp enddo 332 continue @@ -2341,6 +2351,7 @@ c1out open(iout,file=outname,status='unknown') #else if (me.eq.king .or. .not.out1file) & open(iout,file=outname,status='unknown') +#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then write (liczba,'(i3.3)') myrank/nfgtasks @@ -2349,6 +2360,7 @@ c1out open(iout,file=outname,status='unknown') & status='unknown') endif #endif +#undef DEBUG if(me.eq.king) then open(igeom,file=intname,status='unknown',access='append') open(ipdb,file=pdbname,status='unknown') @@ -2734,10 +2746,10 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo @@ -2781,10 +2793,10 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + & write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ", + write (iout,'(a,3i6,f8.2,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif enddo @@ -2811,13 +2823,13 @@ c & fordepth_peak(nhpb_peak+1),npeak ipeak(2,npeak)=i #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak) #else - write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ", + write (iout,'(a,5i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak), & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak), & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak), @@ -2840,11 +2852,11 @@ c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) irestr_type(nhpb)=11 #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + & write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #else - write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,2f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb) #endif @@ -2904,12 +2916,12 @@ c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1) endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + & write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) #else - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) @@ -2938,10 +2950,10 @@ C print *,"in else" endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + write (iout,'(a,4i6,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #endif endif @@ -2984,13 +2996,13 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) dhpb(nhpb)=dist(i+nres*ii,j+nres*jj) #ifdef MPI if (.not.out1file .or. me.eq.king) then - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) endif #else - write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ", + write (iout,'(a,4i6,2f8.2,3f10.5,i5)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb), & irestr_type(nhpb) @@ -3054,8 +3066,8 @@ c & sigma_odl_temp(maxres,maxres,max_template) character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard - integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec, - & ik,iistart,nres_temp + integer ki,i,ii,j,k,l,ii_in_use(maxdim_cont),i_tmp,idomain_tmp, + & irec,ik,iistart,nres_temp integer ilen external ilen logical liiflag,lfirst @@ -3108,7 +3120,7 @@ c endif call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then - write(iout,*) "iset homology_weight " +c write(iout,*) "iset homology_weight " do i=1,homol_nset write(iout,*) i,waga_homology(i) enddo @@ -3122,30 +3134,30 @@ c endif cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) + if (read_homol_frag) then + call read_klapaucjusz + else lim_odl=0 lim_dih=0 c c write(iout,*) 'nnt=',nnt,'nct=',nct c - do i = nnt,nct - do k=1,constr_homology - idomain(k,i)=0 - enddo - enddo - - ii=0 - do i = nnt,nct-2 - do j=i+2,nct - ii=ii+1 - ii_in_use(ii)=0 - enddo - enddo - - if (read_homol_frag) then - call read_klapaucjusz - else +c do i = nnt,nct +c do k=1,constr_homology +c idomain(k,i)=0 +c enddo +c enddo + idomain=0 +c ii=0 +c do i = nnt,nct-2 +c do j=i+2,nct +c ii=ii+1 +c ii_in_use(ii)=0 +c enddo +c enddo + ii_in_use=0 do k=1,constr_homology read(inp,'(a)') pdbfile @@ -3265,8 +3277,8 @@ c & constr_homology endif sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else - ii=ii+1 - l_homo(k,ii)=.false. +c ii=ii+1 +c l_homo(k,ii)=.false. endif enddo enddo @@ -3456,26 +3468,26 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl - write(iout,'(3i5,100(2f8.2,1x,l1,4x))') + write(iout,'(3i7,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii), & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,nct - write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" do i=nnt+2,nct - write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i)), & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" do i=nnt,nct - write(iout,'(i5,100(4f8.2,4x))') i, + write(iout,'(i7,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) enddo @@ -3641,14 +3653,17 @@ c---------------------------------------------------------------------- character*2 kic2 character*24 model_ki_dist, model_ki_angle character*500 controlcard - integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp, - & ik,ll,ii,kk,iistart,iishift,lim_xx + integer ki, i, j, jj,k, l, ii_in_use(maxdim_cont),i_tmp, + & idomain_tmp, + & ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr, + & i01,i10,nnt_chain,nct_chain + integer itype_temp(maxres) double precision distal logical lprn /.true./ integer nres_temp integer ilen external ilen - logical liiflag + logical liiflag,lfirst c c double precision rescore_tmp,x12,y12,z12,rescore2_tmp @@ -3662,14 +3677,48 @@ c For new homol impl c include 'COMMON.VAR' c +c write (iout,*) "READ_KLAPAUCJUSZ" +c print *,"READ_KLAPAUCJUSZ" +c call flush(iout) call getenv("FRAGFILE",fragfile) + write (iout,*) "Opening", fragfile + call flush(iout) open(ientin,file=fragfile,status="old",err=10) - read(ientin,*) constr_homology,nclust - nmodel_start=constr_homology - l_homo = .false. +c write (iout,*) " opened" +c call flush(iout) + sigma_theta=0.0 sigma_d=0.0 sigma_dih=0.0 + l_homo = .false. + + nres_temp=nres + itype_temp=itype + ii=0 + lim_odl=0 + +c write (iout,*) "Entering loop" +c call flush(iout) + + DO IGR = 1,NCHAIN_GROUP + +c write (iout,*) "igr",igr + call flush(iout) + read(ientin,*) constr_homology,nclust + + if (start_from_model) then + nmodel_start=constr_homology + else + nmodel_start=0 + endif + + ii_old=lim_odl + + ichain=iequiv(1,igr) + nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1 + nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1 +c write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain + c Read pdb files do k=1,constr_homology read(ientin,'(a)') pdbfile @@ -3683,10 +3732,10 @@ c Read pdb files stop 34 continue unres_pdb=.false. - nres_temp=nres +c nres_temp=nres call readpdb_template(k) nres_chomo(k)=nres - nres=nres_temp +c nres=nres_temp do i=1,nres rescore(k,i)=0.2d0 rescore2(k,i)=1.0d0 @@ -3705,6 +3754,7 @@ c do ll = 1,ninclust(l) k = inclust(ll,l) +c write (iout,*) "l",l," ll",ll," k",k do i=1,nres idomain(k,i)=0 enddo @@ -3720,7 +3770,7 @@ c Distance restraints c c ... --> odl(k,ii) C Copy the coordinates from reference coordinates (?) - nres_temp=nres +c nres_temp=nres nres=nres_chomo(k) do i=1,2*nres do j=1,3 @@ -3734,11 +3784,13 @@ c write (iout,*) "c(",j,i,") =",c(j,i) thetaref(i)=theta(i) phiref(i)=phi(i) enddo - nres=nres_temp +c nres=nres_temp if (waga_dist.ne.0.0d0) then - ii=0 - do i = nnt,nct-2 - do j=i+2,nct + ii=ii_old +c do i = nnt,nct-2 + do i = nnt_chain,nct_chain-2 +c do j=i+2,nct + do j=i+2,nct_chain x12=c(1,i)-c(1,j) y12=c(2,i)-c(2,j) @@ -3755,9 +3807,9 @@ c write (iout,*) k,i,j,distal,dist2_cut c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", -c & constr_homology - ires_homo(ii)=i - jres_homo(ii)=j +c & constr_homology + ires_homo(ii)=i+chain_border1(1,igr)-1 + jres_homo(ii)=j+chain_border1(1,igr)-1 odl(k,ii)=distal if (read2sigma) then sigma_odl(k,ii)=0 @@ -3793,50 +3845,55 @@ c c Theta, dihedral and SC retraints c if (waga_angle.gt.0.0d0) then - do i = nnt+3,nct + do i = nnt_chain+3,nct_chain + iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then c sigma_dih(k,i)=0.0 cycle endif - dih(k,i)=phiref(i) - sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + dih(k,iii)=phiref(i) + sigma_dih(k,iii)= + & (rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), c & " sigma_dihed",sigma_dih(k,i) - if (sigma_dih(k,i).ne.0) - & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + if (sigma_dih(k,iii).ne.0) + & sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii)) enddo - lim_dih=nct-nnt-2 +c lim_dih=nct-nnt-2 endif if (waga_theta.gt.0.0d0) then - do i = nnt+2,nct + do i = nnt_chain+2,nct_chain + iii=i+chain_border1(1,igr)-1 if (idomain(k,i).eq.0) then c sigma_theta(k,i)=0.0 cycle endif - thetatpl(k,i)=thetaref(i) - sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + thetatpl(k,iii)=thetaref(i) + sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 - if (sigma_theta(k,i).ne.0) - & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,iii).ne.0) + & sigma_theta(k,iii)=1.0d0/ + & (sigma_theta(k,iii)*sigma_theta(k,iii)) enddo endif if (waga_d.gt.0.0d0) then - do i = nnt,nct + do i = nnt_chain,nct_chain + iii=i+chain_border1(1,igr)-1 if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then c sigma_d(k,i)=0.0 cycle endif - xxtpl(k,i)=xxref(i) - yytpl(k,i)=yyref(i) - zztpl(k,i)=zzref(i) - sigma_d(k,i)=rescore(k,i) - if (sigma_d(k,i).ne.0) - & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) - if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 + xxtpl(k,iii)=xxref(i) + yytpl(k,iii)=yyref(i) + zztpl(k,iii)=zzref(i) + sigma_d(k,iii)=rescore(k,i) + if (sigma_d(k,iii).ne.0) + & sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii)) +c if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 enddo endif enddo ! l @@ -3845,41 +3902,101 @@ c c remove distance restraints not used in any model from the list c shift data in all arrays c +c write (iout,*) "ii_old",ii_old if (waga_dist.ne.0.0d0) then - ii=0 +#ifdef DEBUG + write (iout,*) "Distance restraints from templates" + do iii=1,lim_odl + write(iout,'(4i5,100(2f8.2,1x,l1,4x))') + & iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), + & (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), + & ki=1,constr_homology) + enddo +#endif + ii=ii_old liiflag=.true. - do i=nnt,nct-2 - do j=i+2,nct + lfirst=.true. + do i=nnt_chain,nct_chain-2 + do j=i+2,nct_chain ii=ii+1 - if (ii_in_use(ii).eq.0.and.liiflag) then +c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 +c & .and. distal.le.dist2_cut ) then +c write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii) +c call flush(iout) + if (ii_in_use(ii).eq.0.and.liiflag.or. + & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then liiflag=.false. - iistart=ii + i10=ii + if (lfirst) then + lfirst=.false. + iistart=ii + else + if(i10.eq.lim_odl) i10=i10+1 + do ki=0,i10-i01-1 + ires_homo(iistart+ki)=ires_homo(ki+i01) + jres_homo(iistart+ki)=jres_homo(ki+i01) + ii_in_use(iistart+ki)=ii_in_use(ki+i01) + do k=1,constr_homology + odl(k,iistart+ki)=odl(k,ki+i01) + sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01) + l_homo(k,iistart+ki)=l_homo(k,ki+i01) + enddo + enddo + iistart=iistart+i10-i01 + endif endif - if (ii_in_use(ii).ne.0.and..not.liiflag.or. - & .not.liiflag.and.ii.eq.lim_odl) then - if (ii.eq.lim_odl) then - iishift=ii-iistart+1 - else - iishift=ii-iistart - endif + if (ii_in_use(ii).ne.0.and..not.liiflag) then + i01=ii liiflag=.true. - do ki=iistart,lim_odl-iishift - ires_homo(ki)=ires_homo(ki+iishift) - jres_homo(ki)=jres_homo(ki+iishift) - ii_in_use(ki)=ii_in_use(ki+iishift) - do k=1,constr_homology - odl(k,ki)=odl(k,ki+iishift) - sigma_odl(k,ki)=sigma_odl(k,ki+iishift) - l_homo(k,ki)=l_homo(k,ki+iishift) - enddo - enddo - ii=ii-iishift - lim_odl=lim_odl-iishift endif enddo enddo + lim_odl=iistart-1 endif + lll=lim_odl-ii_old + + do i=2,nequiv(igr) + + ichain=iequiv(i,igr) + + do j=nnt_chain,nct_chain + jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr)) + do k=1,constr_homology + dih(k,jj)=dih(k,j) + sigma_dih(k,jj)=sigma_dih(k,j) + thetatpl(k,jj)=thetatpl(k,j) + sigma_theta(k,jj)=sigma_theta(k,j) + xxtpl(k,jj)=xxtpl(k,j) + yytpl(k,jj)=yytpl(k,j) + zztpl(k,jj)=zztpl(k,j) + sigma_d(k,jj)=sigma_d(k,j) + enddo + enddo + + jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr)) +c write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj + do j=ii_old+1,lim_odl + ires_homo(j+lll)=ires_homo(j)+jj + jres_homo(j+lll)=jres_homo(j)+jj + do k=1,constr_homology + odl(k,j+lll)=odl(k,j) + sigma_odl(k,j+lll)=sigma_odl(k,j) + l_homo(k,j+lll)=l_homo(k,j) + enddo + enddo + + ii_old=ii_old+lll + lim_odl=lim_odl+lll + + enddo + + ENDDO ! IGR + + if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2 + nres=nres_temp + itype=itype_temp + return 10 stop "Error in fragment file" end diff --git a/source/unres/src-HCD-5D/refsys.f b/source/unres/src-HCD-5D/refsys.f index 4b7b763..5573f05 100644 --- a/source/unres/src-HCD-5D/refsys.f +++ b/source/unres/src-HCD-5D/refsys.f @@ -21,12 +21,12 @@ c print *,'just initialize' c print *,fail s1=0.0 s2=0.0 - print *,s1,s2 +c print *,s1,s2 do 1 i=1,3 - print *, i2,i3,i4 +c print *, i2,i3,i4 zi=c(i,i2)-c(i,i3) ui=c(i,i4)-c(i,i3) - print *,zi,ui +c print *,zi,ui s1=s1+zi*zi s2=s2+ui*ui z(i)=zi @@ -41,7 +41,7 @@ c print *,fail write(iout,1000) i3,i4,i1 fail=.true. return - print *,'two if pass' +c print *,'two if pass' 4 s1=1.0/s1 s2=1.0/s2 v1=z(2)*u(3)-z(3)*u(2) diff --git a/source/unres/src-HCD-5D/stochfric.F b/source/unres/src-HCD-5D/stochfric.F index c83e9ce..946fb58 100644 --- a/source/unres/src-HCD-5D/stochfric.F +++ b/source/unres/src-HCD-5D/stochfric.F @@ -1,5 +1,8 @@ subroutine friction_force implicit none +#ifdef MPI + include 'mpif.h' +#endif include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' @@ -37,8 +40,12 @@ logical lprn /.false./, checkmode /.false./ #ifdef FIVEDIAG +#ifdef TIMING + include 'COMMON.TIME1' + double precision time01 +#endif c Here accelerations due to friction forces are computed right after forces. - d_t_work=0.0d0 + d_t_work(:6*nres)=0.0d0 do j=1,3 v_work(j,1)=d_t(j,0) v_work(j,nnt)=d_t(j,0) @@ -86,8 +93,14 @@ c inct=chain_border(2,1) write (iout,'(f10.5)') (vvec(i),i=iposc,ind) #endif c write (iout,*) "chain",i," ind",ind," n",n +#ifdef TIMING + time01=MPI_Wtime() +#endif call fivediagmult(n,DMfric(iposc),DU1fric(iposc), & DU2fric(iposc),vvec(iposc),rs) +#ifdef TIMING + time_fricmatmult=time_fricmatmult+MPI_Wtime()-time01 +#endif #ifdef DEBUG write (iout,*) "rs" write (iout,'(f10.5)') (rs(i),i=1,n) @@ -104,7 +117,7 @@ c & "index",3*(i-1)+j,"rs",rs(i-iposc+1) write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3) #endif #else - do i=0,MAXRES2 + do i=0,2*nres do j=1,3 friction(j,i)=0.0d0 enddo @@ -277,7 +290,7 @@ c----------------------------------------------------- integer ichain,innt,inct,iposc #endif - do i=0,MAXRES2 + do i=0,2*nres do j=1,3 stochforc(j,i)=0.0d0 enddo @@ -567,9 +580,9 @@ C gamsc(ntyp1)=1.0d0 enddo endif #ifdef FIVEDIAG - DMfric=0.0d0 - DU1fric=0.0d0 - DU2fric=0.0d0 + DMfric(:2*nres)=0.0d0 + DU1fric(:2*nres)=0.0d0 + DU2fric(:2*nres)=0.0d0 ind=1 do ichain=1,nchain innt=chain_border(1,ichain) diff --git a/source/unres/src-HCD-5D/timing.F b/source/unres/src-HCD-5D/timing.F index 7bd51b8..98ce59a 100644 --- a/source/unres/src-HCD-5D/timing.F +++ b/source/unres/src-HCD-5D/timing.F @@ -50,6 +50,18 @@ C Calculate the initial time, if it is not zero (e.g. for the SUN). time_fricmatmult=0.0d0 time_fsample=0.0d0 time_SAXS=0.0d0 + time_list=0.0d0 + time_evdw=0.0d0 + time_evdw_short=0.0d0 + time_evdw_long=0.0d0 + time_eelec=0.0d0 + time_eelec_short=0.0d0 + time_eelec_long=0.0d0 + time_escp=0.0d0 + time_escp_short=0.0d0 + time_escp_long=0.0d0 + time_escpsetup=0.0d0 + time_escpcalc=0.0d0 #endif cd print *,' in SET_TIMERS stime=',stime return @@ -287,6 +299,7 @@ C--------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.TIME1' include 'COMMON.SETUP' + include 'COMMON.MD' #ifdef MPI time1=MPI_WTIME() write (iout,'(80(1h=)/a/(80(1h=)))') @@ -318,6 +331,30 @@ C--------------------------------------------------------------------------- & time_bcast+time_reduce+time_gather+time_scatter+ & time_sendrecv+time_barrier_g+time_barrier_e+time_bcastc write (*,*) "Processor",fg_rank,myrank," enecalc",time_enecalc +#ifdef TIMING_ENE + if (RESPA) then + write (*,*) "Processor",fg_rank,myrank," evdw_long", + & time_evdw_long + write (*,*) "Processor",fg_rank,myrank," evdw_short", + & time_evdw_short + write (*,*) "Processor",fg_rank,myrank," eelec_long", + & time_eelec_long + write (*,*) "Processor",fg_rank,myrank," eelec_short", + & time_eelec_short + write (*,*) "Processor",fg_rank,myrank," escp_long", + & time_escp_long + write (*,*) "Processor",fg_rank,myrank," escp_short", + & time_escp_short + else + write (*,*) "Processor",fg_rank,myrank," evdw",time_evdw + write (*,*) "Processor",fg_rank,myrank," eelec",time_eelec + write (*,*) "Processor",fg_rank,myrank," escp",time_escp + write (*,*) "Processor",fg_rank,myrank," escpsetup", + & time_escpsetup + write (*,*) "Processor",fg_rank,myrank," escpcalc", + & time_escpcalc + endif +#endif write (*,*) "Processor",fg_rank,myrank," sumene",time_sumene write (*,*) "Processor",fg_rank,myrank," intfromcart", & time_intfcart diff --git a/source/unres/src-HCD-5D/unres.F b/source/unres/src-HCD-5D/unres.F index 978fa59..2e0ebaf 100644 --- a/source/unres/src-HCD-5D/unres.F +++ b/source/unres/src-HCD-5D/unres.F @@ -213,11 +213,11 @@ c include 'COMMON.CONTACTS' include 'COMMON.REMD' include 'COMMON.MD' include 'COMMON.SBRIDGE' - integer i,icall,iretcode,nfun + integer i,it,icall,iretcode,nfun common /srutu/ icall integer nharp,iharp(4,maxres/3) integer nft_sc - logical fail + logical fail,secondary_str /.true./ double precision energy(0:n_ene),etot,etota double precision energy_long(0:n_ene),energy_short(0:n_ene) double precision rms,frac,frac_nn,co @@ -243,6 +243,24 @@ c include 'COMMON.CONTACTS' write (iout,*) "Energy evaluation/minimization" call chainbuild_cart c print *,'dc',dc(1,0),dc(2,0),dc(3,0) + if (nran_start.gt.0) then + write (iout,*) + & "Chains will be regenerated starting from residue",nran_start + do it=1,100 + call gen_rand_conf_mchain(nran_start,*10) + write (iout,*) "Conformation successfully generated",it + goto 11 + 10 write (iout,*) "Problems with regenerating chains",it + enddo + 11 continue + write (iout,*) "Cartesian coords after chain rebuild" + call cartprint + call chainbuild_cart + write (iout,*) "Cartesian coords after chainbuild_ecart" + call cartprint + call int_from_cart1(.false.) + call intout + endif if (split_ene) then print *,"Processor",myrank," after chainbuild" icall=1 @@ -274,9 +292,11 @@ c print *,"after etotal" etota = energy(0) etot =etota call enerprint(energy(0)) + if (secondary_str) then call hairpin(.true.,nharp,iharp) c print *,'after hairpin' call secondary2(.true.) + endif c print *,'after secondary' if (minim) then crc overlap test @@ -292,6 +312,9 @@ crc overlap test write (iout,*) 'Calling OVERLAP_SC' call overlap_sc(fail) write (iout,*) "After overlap_sc" +c cartname=prefix(:ilen(prefix))//'.x' +c potE=etot +c call cartoutx(0.0d0) endif if (searchsc) then @@ -330,10 +353,12 @@ c print *,'Calling MINIMIZE.' #endif print *,'# eval/s',evals print *,'refstr=',refstr + if (secondary_str) then call hairpin(.false.,nharp,iharp) - print *,'after hairpin' +c print *,'after hairpin' call secondary2(.true.) - print *,'after secondary' +c print *,'after secondary' + endif call etotal(energy(0)) etot = energy(0) call enerprint(energy(0)) -- 1.7.9.5