From a937b3a5ce9b4ee672e131cff3e8d022ca237b9e Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Thu, 26 Mar 2020 18:49:14 +0100 Subject: [PATCH] Adam's update --- source/unres/src-HCD-5D/COMMON.INTERACT | 10 + source/unres/src-HCD-5D/COMMON.MD | 4 +- source/unres/src-HCD-5D/MD_A-MTS.F | 30 +- .../unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos | 4 +- source/unres/src-HCD-5D/energy_p_new_barrier.F | 132 ++++-- .../unres/src-HCD-5D/energy_p_new_barrier.F.safe | 221 +++++---- source/unres/src-HCD-5D/initialize_p.F | 3 +- source/unres/src-HCD-5D/lagrangian_lesyng.F | 14 +- source/unres/src-HCD-5D/make_xx_list.F | 489 ++++++++++++++++++++ source/unres/src-HCD-5D/readrtns_CSA.F | 5 + source/unres/src-HCD-5D/stochfric.F | 7 + 11 files changed, 778 insertions(+), 141 deletions(-) create mode 100644 source/unres/src-HCD-5D/make_xx_list.F diff --git a/source/unres/src-HCD-5D/COMMON.INTERACT b/source/unres/src-HCD-5D/COMMON.INTERACT index 14d92ef..3440239 100644 --- a/source/unres/src-HCD-5D/COMMON.INTERACT +++ b/source/unres/src-HCD-5D/COMMON.INTERACT @@ -21,6 +21,16 @@ & iscpstart(maxres,maxint_gr),iscpend(maxres,maxint_gr), & 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 + integer newcontlisti(200*maxres),newcontlistj(200*maxres), + & newcontlistppi(200*maxres),newcontlistppj(200*maxres), + & newcontlistscpi(200*maxres),newcontlistscpj(200*maxres), + & g_listscsc_start,g_listscsc_end,g_listpp_start,g_listpp_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,newcontlistscpi,newcontlistscpj,g_listscp_start, + & g_listscp_end 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 6988bd8..cea18eb 100644 --- a/source/unres/src-HCD-5D/COMMON.MD +++ b/source/unres/src-HCD-5D/COMMON.MD @@ -3,12 +3,12 @@ & dvmax,damax,edriftmax integer n_timestep,ntwx,ntwe,lang,count_reset_moment, & count_reset_vel,ntime_split,ntime_split0, - & maxtime_split + & maxtime_split,itime_mat,imatupdate 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, + & ntime_split,ntime_split0,maxtime_split,itime_mat,imatupdate, & ntwx,ntwe,lang,large,print_compon,tbf,rest,preminim, & reset_moment,reset_vel,count_reset_moment,count_reset_vel, & rattle,RESPA diff --git a/source/unres/src-HCD-5D/MD_A-MTS.F b/source/unres/src-HCD-5D/MD_A-MTS.F index 08852ba..4498b4a 100644 --- a/source/unres/src-HCD-5D/MD_A-MTS.F +++ b/source/unres/src-HCD-5D/MD_A-MTS.F @@ -166,14 +166,14 @@ c Entering the MD loop & .and. mod(itime,count_reset_vel).eq.0) then call random_vel write(iout,'(a,f20.2)') - & "Velocities reset to random values, time",totT + & "Velocities reset to random values, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=d_t(j,i) enddo enddo endif - if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then + if (reset_moment .and. mod(itime,count_reset_moment).eq.0) then call inertia_tensor call vcm_vel(vcm) do j=1,3 @@ -182,7 +182,7 @@ c Entering the MD loop call kinetic(EK) kinetic_T=2.0d0/(dimen3*Rb)*EK scalfac=dsqrt(T_bath/kinetic_T) - write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT + write(iout,'(a,f20.2)') "Momenta zeroed out, time",totT do i=0,2*nres do j=1,3 d_t_old(j,i)=scalfac*d_t(j,i) @@ -209,6 +209,7 @@ c Variable time step algorithm. stop #endif endif + itime_mat=itime if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0) then call statout(itime) @@ -1811,10 +1812,11 @@ c rest2name = prefix(:ilen(prefix))//'.rst' write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), & (d_t(j,i+nres),j=1,3) enddo + endif + if (lang.eq.0) then c Zeroing the total angular momentum of the system write(iout,*) "Calling the zero-angular & momentum subroutine" - endif call inertia_tensor c Getting the potential energy and forces and velocities and accelerations call vcm_vel(vcm) @@ -1829,8 +1831,14 @@ c Removing the velocity of the center of mass if(me.eq.king.or..not.out1file)then write (iout,*) "vcm right after adjustment:" write (iout,*) (vcm(j),j=1,3) + write (iout,*) "Initial velocities after adjustment" + do i=0,nres + write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), + & (d_t(j,i+nres),j=1,3) + enddo + call flush(iout) + endif endif - call flush(iout) write (iout,*) "init_MD before initial structure REST ",rest if (.not.rest) then 122 continue @@ -2050,6 +2058,7 @@ c write(iout,*) (potEcomp(i),i=0,n_ene) c write (iout,*) "PotE-homology",potE totE=EK+potE itime=0 + itime_mat=itime if (ntwe.ne.0) call statout(itime) if(me.eq.king.or..not.out1file) & write (iout,'(/a/3(a25,1pe14.5/))') "Initial:", @@ -2351,6 +2360,9 @@ c write (iout,*) "ichain",ichain," innt",innt," inct",inct ! #define WLOS #ifdef WLOS + if (nnt.eq.1) then + d_t(:,0)=d_t(:,1) + endif do i=1,nres if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then do j=1,3 @@ -2363,11 +2375,17 @@ c write (iout,*) "ichain",ichain," innt",innt," inct",inct enddo end if enddo + d_t(:,nres)=0.0d0 d_t(:,nct)=0.0d0 + d_t(:,2*nres)=0.0d0 + if (nnt.gt.1) then + d_t(:,0)=d_t(:,1) + d_t(:,1)=0.0d0 + endif c d_a(:,0)=d_a(:,1) c d_a(:,1)=0.0d0 c write (iout,*) "Shifting accelerations" - do ichain=1,nchain + do ichain=2,nchain c write (iout,*) "ichain",chain_border1(1,ichain)-1, c & chain_border1(1,ichain) d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain)) diff --git a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos index 5485424..32a0dec 100644 --- a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos +++ b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos @@ -38,7 +38,7 @@ object = unres.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \ matmult.o readrtns_CSA.o parmread.o gen_rand_conf.o printmat.o map.o \ pinorm.o randgens.o rescode.o intcor.o timing.o misc.o \ cart2intgrad.o checkder_p.o contact_cp econstr_local.o econstr_qlike.o \ - econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o \ + econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o make_xx_list \ energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \ cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \ mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o \ @@ -95,7 +95,7 @@ NEWCORR: ${object} xdrf/libxdrf.a NEWCORR5D: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \ -DSPLITELE -DLANG0 -DNEWCORR -DCORRCD -DFIVEDIAG -DLBFGS #-DMYGAUSS #-DTIMING -NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5.exe +NEWCORR5D: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD5-l.exe NEWCORR5D: ${object_lbfgs} ${object} fdisy.o fdiag.o machpd.o kinetic_CASC.o xdrf/libxdrf.a gcc -o compinfo compinfo.c ./compinfo | true 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 2a588bd..ef19809 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -30,6 +30,7 @@ c include 'COMMON.MD' include 'COMMON.SPLITELE' include 'COMMON.TORCNSTR' include 'COMMON.SAXS' + include 'COMMON.MD' 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, @@ -117,6 +118,12 @@ c call chainbuild_cart edfanei=0.0d0 edfabet=0.0d0 #endif + if (nfgtasks.gt.1) then + call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR) + endif + if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list + if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list + if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list c print *,'Processor',myrank,' calling etotal ipot=',ipot c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct #else @@ -1460,14 +1467,17 @@ C #endif double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,num_conti,iint + integer i,j,k,itypi,itypj,itypi1,num_conti,iint,icont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision fcont,fprimcont double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1479,10 +1489,10 @@ C Change 12/1/95 C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),iend(i,iint) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi @@ -1587,8 +1597,8 @@ cd & i,j,(gacont(kk,num_conti,i),kk=1,3) endif endif #endif - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint C Change 12/1/95 #ifdef FOURBODY num_cont(i)=num_conti @@ -1630,14 +1640,17 @@ C include 'COMMON.SPLITELE' double precision gg(3) double precision evdw,evdwij - integer i,j,k,itypi,itypj,itypi1,iint + integer i,j,k,itypi,itypj,itypi1,iint,icont double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck double precision sscale,sscagrad c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1647,8 +1660,8 @@ c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi @@ -1695,8 +1708,8 @@ cgrad do l=1,3 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i do i=1,nct do j=1,3 @@ -1727,7 +1740,7 @@ C integer icall common /srutu/ icall double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,icont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, & sss1,sssgrad1 double precision sscale,sscagrad @@ -1742,7 +1755,10 @@ c else lprn=.false. c endif ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1757,8 +1773,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -1835,8 +1851,8 @@ C Calculate radial part of the gradient C Calculate the angular part of the gradient and sum add the contributions C to the appropriate components of the Cartesian gradient. call sc_grad - enddo ! j - enddo ! iint +! enddo ! j +! enddo ! iint enddo ! i c stop return @@ -1864,7 +1880,7 @@ C logical lprn integer xshift,yshift,zshift,subchap double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,icont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, @@ -1882,7 +1898,10 @@ C we have the original box) C do xshift=-1,1 C do yshift=-1,1 C do zshift=-1,1 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -1963,8 +1982,8 @@ c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN c write(iout,*) "PRZED ZWYKLE", evdwij @@ -2211,8 +2230,8 @@ C Calculate angular part of the gradient. c call sc_grad_scale(sss) call sc_grad ENDIF ! dyn_ss - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i C enddo ! zshift C enddo ! yshift @@ -2244,7 +2263,7 @@ C common /srutu/ icall logical lprn double precision evdw - integer itypi,itypj,itypi1,iint,ind + integer itypi,itypj,itypi1,iint,ind,icont double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij, & xi,yi,zi,fac_augm,e_augm double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, @@ -2257,7 +2276,10 @@ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon lprn=.false. c if (icall.eq.0) lprn=.true. ind=0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -2307,8 +2329,8 @@ c dsci_inv=dsc_inv(itypi) C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) - do j=istart(i,iint),iend(i,iint) +c do iint=1,nint_gr(i) +c do j=istart(i,iint),iend(i,iint) ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -2460,8 +2482,8 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. c call sc_grad_scale(sss) call sc_grad - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i end C----------------------------------------------------------------------------- @@ -2612,7 +2634,10 @@ c include 'COMMON.CONTACTS' dimension gg(3) cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct evdw=0.0D0 - do i=iatsc_s,iatsc_e +c do i=iatsc_s,iatsc_e + do icont=g_listscsc_start,g_listscsc_end + i=newcontlisti(icont) + j=newcontlistj(icont) itypi=iabs(itype(i)) if (itypi.eq.ntyp1) cycle itypi1=iabs(itype(i+1)) @@ -2622,10 +2647,10 @@ cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct C C Calculate SC interaction energy. C - do iint=1,nint_gr(i) +c do iint=1,nint_gr(i) cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), cd & 'iend=',iend(i,iint) - do j=istart(i,iint),iend(i,iint) +c do j=istart(i,iint),iend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi @@ -2661,8 +2686,8 @@ cgrad do l=1,3 cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) cgrad enddo cgrad enddo - enddo ! j - enddo ! iint +c enddo ! j +c enddo ! iint enddo ! i return end @@ -3843,7 +3868,10 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c CTU KURWA - do i=iatel_s,iatel_e +c do i=iatel_s,iatel_e + do icont=g_listpp_start,g_listpp_end + i=newcontlistppi(icont) + j=newcontlistppj(icont) C do i=75,75 c if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 @@ -3904,7 +3932,7 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) #endif C I TU KURWA - do j=ielstart(i),ielend(i) +c do j=ielstart(i),ielend(i) C do j=16,17 C write (iout,*) i,j C if (j.le.1) cycle @@ -3917,7 +3945,7 @@ c & .or.itype(j+2).eq.ntyp1 c & .or.itype(j-1).eq.ntyp1 &) cycle call eelecij(i,j,ees,evdw1,eel_loc) - enddo ! j +c enddo ! j #ifdef FOURBODY num_cont_hb(i)=num_conti #endif @@ -5546,7 +5574,10 @@ cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e C do xshift=-1,1 C do yshift=-1,1 C do zshift=-1,1 - do i=iatscp_s,iatscp_e +c do i=iatscp_s,iatscp_e + do icont=g_listscp_start,g_listscp_end + i=newcontlistscpi(icont) + j=newcontlistscpj(icont) 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)) @@ -5586,9 +5617,9 @@ c endif C xi=xi+xshift*boxxsize C yi=yi+yshift*boxysize C zi=zi+zshift*boxzsize - do iint=1,nscp_gr(i) +c do iint=1,nscp_gr(i) - do j=iscpstart(i,iint),iscpend(i,iint) +c do j=iscpstart(i,iint),iscpend(i,iint) if (itype(j).eq.ntyp1) cycle itypj=iabs(itype(j)) C Uncomment following three lines for SC-p interactions @@ -5711,9 +5742,9 @@ cgrad enddo gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo - enddo +c enddo - enddo ! iint +c enddo ! iint enddo ! i C enddo !zshift C enddo !yshift @@ -5741,7 +5772,7 @@ C include 'COMMON.SPLITELE' integer xshift,yshift,zshift double precision ggg(3) - integer i,iint,j,k,iteli,itypj,subchap + integer i,iint,j,k,iteli,itypj,subchap,icont double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, & fac,e1,e2,rij double precision evdw2,evdw2_14,evdwij @@ -5757,7 +5788,10 @@ C do xshift=-1,1 C do yshift=-1,1 C do zshift=-1,1 if (energy_dec) write (iout,*) "escp:",r_cut_int,rlamb - do i=iatscp_s,iatscp_e +c do i=iatscp_s,iatscp_e + do icont=g_listscp_start,g_listscp_end + i=newcontlistscpi(icont) + j=newcontlistscpj(icont) 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)) @@ -5800,9 +5834,9 @@ c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or. c & (zi.lt.((zshift-0.5d0)*boxzsize))) then c go to 136 c endif - do iint=1,nscp_gr(i) +c do iint=1,nscp_gr(i) - do j=iscpstart(i,iint),iscpend(i,iint) +c do j=iscpstart(i,iint),iscpend(i,iint) itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle C Uncomment following three lines for SC-p interactions @@ -5937,9 +5971,9 @@ cgrad enddo gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) enddo c endif !endif for sscale cutoff - enddo ! j +c enddo ! j - enddo ! iint +c enddo ! iint enddo ! i c enddo !zshift c enddo !yshift diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F.safe b/source/unres/src-HCD-5D/energy_p_new_barrier.F.safe index ae8e449..2a588bd 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F.safe +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F.safe @@ -66,7 +66,8 @@ C FG slaves as WEIGHTS array. weights_(17)=wbond weights_(18)=scal14 weights_(21)=wsccor - weights_(22)=wtube + weights_(22)=wliptran + weights_(25)=wtube weights_(26)=wsaxs weights_(28)=wdfa_dist weights_(29)=wdfa_tor @@ -98,7 +99,8 @@ C FG slaves receive the WEIGHTS array wbond=weights(17) scal14=weights(18) wsccor=weights(21) - wtube=weights(22) + wliptran=weights(22) + wtube=weights(25) wsaxs=weights(26) wdfa_dist=weights_(28) wdfa_tor=weights_(29) @@ -387,6 +389,8 @@ C based on partition function C print *,"przed lipidami" if (wliptran.gt.0) then call Eliptransfer(eliptran) + else + eliptran=0.0d0 endif C print *,"za lipidami" if (AFMlog.gt.0) then @@ -1449,6 +1453,7 @@ C include 'COMMON.SBRIDGE' include 'COMMON.NAMES' include 'COMMON.IOUNITS' + include 'COMMON.SPLITELE' #ifdef FOURBODY include 'COMMON.CONTACTS' include 'COMMON.CONTMAT' @@ -1457,8 +1462,9 @@ C double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,num_conti,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, - & sigij,r0ij,rcut + & sigij,r0ij,rcut,sqrij,sss1,sssgrad1 double precision fcont,fprimcont + double precision sscale,sscagrad c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e @@ -1485,6 +1491,11 @@ cd & 'iend=',iend(i,iint) C Change 12/1/95 to calculate four-body interactions rij=xj*xj+yj*yj+zj*zj rrij=1.0D0/rij + sqrij=dsqrt(rij) + sss1=sscale(sqrij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(sqrij,r_cut_int) + c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 @@ -1498,11 +1509,12 @@ cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj), cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - evdw=evdw+evdwij + evdw=evdw+sss1*evdwij C C Calculate the components of the gradient in DC and X C - fac=-rrij*(e1+evdwij) + fac=-rrij*(e1+evdwij)*sss1 + & +evdwij*sssgrad1/sqrij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1615,12 +1627,14 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.NAMES' + include 'COMMON.SPLITELE' double precision gg(3) double precision evdw,evdwij integer i,j,k,itypi,itypj,itypi1,iint double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij, - & fac_augm,e_augm,r_inv_ij,r_shift_inv + & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1 logical scheck + double precision sscale,sscagrad c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e @@ -1645,6 +1659,9 @@ C e_augm=augm(itypi,itypj)*fac_augm r_inv_ij=dsqrt(rrij) rij=1.0D0/r_inv_ij + sss1=sscale(rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(rij,r_cut_int) r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj)) fac=r_shift_inv**expon C have you changed here? @@ -1658,11 +1675,12 @@ cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, cd & (c(k,i),k=1,3),(c(k,j),k=1,3) - evdw=evdw+evdwij + evdw=evdw+evdwij*sss1 C C Calculate the components of the gradient in DC and X C fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2) + & +evdwij*sssgrad1*r_inv_ij/expon gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac @@ -1705,11 +1723,14 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SPLITELE' integer icall common /srutu/ icall double precision evdw integer itypi,itypj,itypi1,iint,ind - double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi + double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi, + & sss1,sssgrad1 + double precision sscale,sscagrad c double precision rrsave(maxdim) logical lprn evdw=0.0D0 @@ -1775,6 +1796,9 @@ cd else cd rrij=rrsave(ind) cd endif rij=dsqrt(rrij) + sss1=sscale(1.0d0/rij,r_cut_int) + if (sss1.eq.0.0d0) cycle + sssgrad1=sscagrad(1.0d0/rij,r_cut_int) C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions @@ -1787,7 +1811,7 @@ C have you changed here? eps2der=evdwij*eps3rt eps3der=evdwij*eps2rt evdwij=evdwij*eps2rt*eps3rt - evdw=evdw+evdwij + evdw=evdw+sss1*evdwij if (lprn) then sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa @@ -1803,6 +1827,7 @@ C Calculate gradient components. fac=-expon*(e1+evdwij) sigder=fac/sigsq fac=rrij*fac + & +evdwij*sssgrad1/sss1*rij C Calculate radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac @@ -2108,11 +2133,10 @@ c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) sss=sscale(1.0d0/rij,r_cut_int) - sssgrad=sscagrad(1.0d0/rij,r_cut_int) - c write (iout,'(a7,4f8.3)') c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb - if (sss.gt.0.0d0) then + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -2159,8 +2183,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 & evdwij endif - if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') - & 'evdw',i,j,evdwij + if (energy_dec) write (iout,'(a,2i5,3f10.5)') + & 'r sss evdw',i,j,rij,sss,evdwij C Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -2169,13 +2193,13 @@ C Calculate gradient components. fac=rij*fac c print '(2i4,6f8.4)',i,j,sss,sssgrad* c & evdwij,fac,sigma(itypi,itypj),expon - fac=fac+evdwij/sss*sssgrad*rij + fac=fac+evdwij*sssgrad/sss*rij c fac=0.0d0 C Calculate the radial part of the gradient gg_lipi(3)=eps1*(eps2rt*eps2rt) - &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* - & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) - &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) + & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip* + & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj)) + & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))) gg_lipj(3)=ssgradlipj*gg_lipi(3) gg_lipi(3)=gg_lipi(3)*ssgradlipi C gg_lipi(3)=0.0d0 @@ -2184,8 +2208,8 @@ C gg_lipj(3)=0.0d0 gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. +c call sc_grad_scale(sss) call sc_grad - endif ENDIF ! dyn_ss enddo ! j enddo ! iint @@ -2214,6 +2238,7 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SPLITELE' integer xshift,yshift,zshift,subchap integer icall common /srutu/ icall @@ -2224,7 +2249,7 @@ C & xi,yi,zi,fac_augm,e_augm double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij, & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe, - & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip + & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1 double precision dist,sscale,sscagrad,sscagradlip,sscalelip evdw=0.0D0 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon @@ -2384,6 +2409,9 @@ C write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj dzj=dc_norm(3,nres+j) rrij=1.0D0/(xj*xj+yj*yj+zj*zj) rij=dsqrt(rrij) + sss=sscale(1.0d0/rij,r_cut_int) + if (sss.eq.0.0d0) cycle + sssgrad=sscagrad(1.0d0/rij,r_cut_int) C Calculate angle-dependent terms of energy and contributions to their C derivatives. call sc_angular @@ -2424,12 +2452,13 @@ C Calculate gradient components. fac=-expon*(e1+evdwij)*rij_shift sigder=fac*sigder fac=rij*fac-2*expon*rrij*e_augm - fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij + fac=fac+(evdwij+e_augm)*sssgrad/sss*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac C Calculate angular part of the gradient. +c call sc_grad_scale(sss) call sc_grad enddo ! j enddo ! iint @@ -3909,7 +3938,7 @@ cd print *,"Processor",fg_rank," t_eelecij",t_eelecij end C------------------------------------------------------------------------------- subroutine eelecij(i,j,ees,evdw1,eel_loc) - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -3933,14 +3962,33 @@ C------------------------------------------------------------------------------- include 'COMMON.TIME1' include 'COMMON.SPLITELE' include 'COMMON.SHIELD' - dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), + double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4), & gmuij2(4),gmuji2(4) + double precision dxi,dyi,dzi + double precision dx_normi,dy_normi,dz_normi,aux + integer j1,j2,lll,num_conti common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, & num_conti,j1,j2 + integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ilist,iresshield + double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp + double precision ees,evdw1,eel_loc,aaa,bbb,ael3i + double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj, + & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4, + & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa, + & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der, + & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij, + & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp, + & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp, + & ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield + double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji + double precision dist_init,xj_safe,yj_safe,zj_safe, + & xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi + double precision sscale,sscagrad,scalar + c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions #ifdef MOMENT double precision scal_el /1.0d0/ @@ -4044,8 +4092,9 @@ C yj=yj-ymedi C zj=zj-zmedi rij=xj*xj+yj*yj+zj*zj - sss=sscale(sqrt(rij),r_cut_int) - sssgrad=sscagrad(sqrt(rij),r_cut_int) + sss=sscale(dsqrt(rij),r_cut_int) + if (sss.eq.0.0d0) return + sssgrad=sscagrad(dsqrt(rij),r_cut_int) c if (sss.gt.0.0d0) then rrmij=1.0D0/rij rij=dsqrt(rij) @@ -4080,7 +4129,7 @@ C fac_shield(j)=0.6 fac_shield(i)=1.0 fac_shield(j)=1.0 eesij=(el1+el2) - ees=ees+eesij + ees=ees+eesij*sss endif evdw1=evdw1+evdwij*sss cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)') @@ -4089,11 +4138,10 @@ cd & 1.0D0/dsqrt(rrmij),evdwij,eesij, cd & xmedi,ymedi,zmedi,xj,yj,zj if (energy_dec) then - write (iout,'(a6,2i5,0pf7.3,2i5,3e11.3)') - &'evdw1',i,j,evdwij - &,iteli,itelj,aaa,evdw1,sss - write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, - &fac_shield(i),fac_shield(j) + write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,3f10.5)') + & 'evdw1',i,j,evdwij,iteli,itelj,aaa,evdw1,sss,rij + write (iout,'(a6,2i5,0pf7.3,2f8.3)') 'ees',i,j,eesij, + & fac_shield(i),fac_shield(j) endif C @@ -4110,9 +4158,10 @@ C * * Radial derivatives. First process both termini of the fragment (i,j) * - ggg(1)=facel*xj - ggg(2)=facel*yj - ggg(3)=facel*zj + aux=facel*sss+rmij*sssgrad*eesij + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then C print *,i,j @@ -4146,10 +4195,10 @@ C endif iresshield=shield_list(ilist,j) do k=1,3 rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j) - & *2.0 + & *2.0*sss gshieldx(k,iresshield)=gshieldx(k,iresshield)+ & rlocshield - & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0 + & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield C & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j) @@ -4172,13 +4221,13 @@ C endif do k=1,3 gshieldc(k,i)=gshieldc(k,i)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j)=gshieldc(k,j)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss gshieldc(k,i-1)=gshieldc(k,i-1)+ - & grad_shield(k,i)*eesij/fac_shield(i)*2.0 + & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss gshieldc(k,j-1)=gshieldc(k,j-1)+ - & grad_shield(k,j)*eesij/fac_shield(j)*2.0 + & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss enddo endif @@ -4209,15 +4258,10 @@ cgrad do l=1,3 cgrad gelc(l,k)=gelc(l,k)+ggg(l) cgrad enddo cgrad enddo - if (sss.gt.0.0) then - ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj - ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj - ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj - else - ggg(1)=0.0 - ggg(2)=0.0 - ggg(3)=0.0 - endif + facvdw=facvdw+sssgrad*rmij*evdwij + ggg(1)=facvdw*xj + ggg(2)=facvdw*yj + ggg(3)=facvdw*zj c do k=1,3 c ghalf=0.5D0*ggg(k) c gvdwpp(k,i)=gvdwpp(k,i)+ghalf @@ -4238,10 +4282,11 @@ cgrad enddo cgrad enddo #else C MARYSIA - facvdw=(ev1+evdwij)*sss + facvdw=(ev1+evdwij) facel=(el1+eesij) fac1=fac - fac=-3*rrmij*(facvdw+facvdw+facel) + fac=-3*rrmij*(facvdw+facvdw+facel)*sss + & +(evdwij+eesij)*sssgrad*rrmij erij(1)=xj*rmij erij(2)=yj*rmij erij(3)=zj*rmij @@ -4297,7 +4342,7 @@ cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3), cd & (dcosg(k),k=1,3) do k=1,3 ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))* - & fac_shield(i)**2*fac_shield(j)**2 + & fac_shield(i)**2*fac_shield(j)**2*sss enddo c do k=1,3 c ghalf=0.5D0*ggg(k) @@ -4317,11 +4362,11 @@ C print *,"before22", gelc_long(1,i), gelc_long(1,j) do k=1,3 gelc(k,i)=gelc(k,i) & +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) - & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)) + & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss & *fac_shield(i)**2*fac_shield(j)**2 gelc(k,j)=gelc(k,j) & +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) - & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)) + & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))*sss & *fac_shield(i)**2*fac_shield(j)**2 gelc_long(k,j)=gelc_long(k,j)+ggg(k) gelc_long(k,i)=gelc_long(k,i)-ggg(k) @@ -4557,7 +4602,7 @@ C fac_shield(i)=0.4 C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') c & 'eelloc',i,j,eel_loc_ij C Now derivative over eel_loc @@ -4615,7 +4660,7 @@ C Calculate patrial derivative for theta angle & +a23*gmuij1(2) & +a32*gmuij1(3) & +a33*gmuij1(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c write(iout,*) "derivative over thatai" c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), c & a33*gmuij1(4) @@ -4631,7 +4676,7 @@ c & a33*gmuij2(4) & +a33*gmuij2(4) gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+ & geel_loc_ij*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss c Derivative over j residue geel_loc_ji=a22*gmuji1(1) @@ -4644,7 +4689,7 @@ c & a33*gmuji1(4) gloc(nphi+j,icg)=gloc(nphi+j,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss geel_loc_ji= & +a22*gmuji2(1) @@ -4656,7 +4701,7 @@ c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), c & a33*gmuji2(4) gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+ & geel_loc_ji*wel_loc - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss #endif cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij @@ -4672,17 +4717,21 @@ C Partial derivatives in virtual-bond dihedral angles gamma & gel_loc_loc(i-1)=gel_loc_loc(i-1)+ & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) + aux=eel_loc_ij/sss*sssgrad*rmij + ggg(1)=aux*xj + ggg(2)=aux*yj + ggg(3)=aux*zj do l=1,3 - ggg(l)=(agg(l,1)*muij(1)+ + ggg(l)=ggg(l)+(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) cgrad ghalf=0.5d0*ggg(l) @@ -4698,19 +4747,19 @@ C Remaining derivatives of eello do l=1,3 gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss enddo ENDIF @@ -4800,9 +4849,9 @@ C fac_shield(i)=0.4d0 C fac_shield(j)=0.6d0 endif ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss C Diagnostics. Comment out or remove after debugging! c ees0p(num_conti,i)=0.5D0*fac3*ees0pij c ees0m(num_conti,i)=0.5D0*fac3*ees0mij @@ -4851,11 +4900,17 @@ cd fprimcont=0.0D0 gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k) enddo gggp(1)=gggp(1)+ees0pijp*xj + & +ees0p(num_conti,i)/sss*rmij*xj*sssgrad gggp(2)=gggp(2)+ees0pijp*yj + & +ees0p(num_conti,i)/sss*rmij*yj*sssgrad gggp(3)=gggp(3)+ees0pijp*zj + & +ees0p(num_conti,i)/sss*rmij*zj*sssgrad gggm(1)=gggm(1)+ees0mijp*xj + & +ees0m(num_conti,i)/sss*rmij*xj*sssgrad gggm(2)=gggm(2)+ees0mijp*yj + & +ees0m(num_conti,i)/sss*rmij*yj*sssgrad gggm(3)=gggm(3)+ees0mijp*zj + & +ees0m(num_conti,i)/sss*rmij*zj*sssgrad C Derivatives due to the contact function gacont_hbr(1,num_conti,i)=fprimcont*xj gacont_hbr(2,num_conti,i)=fprimcont*yj @@ -4870,28 +4925,28 @@ cgrad ghalfm=0.5D0*gggm(k) gacontp_hb1(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb2(k,num_conti,i)=!ghalfp & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontp_hb3(k,num_conti,i)=gggp(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb1(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i)) & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb2(k,num_conti,i)=!ghalfm & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j)) & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss gacontm_hb3(k,num_conti,i)=gggm(k) - & *fac_shield(i)*fac_shield(j) + & *fac_shield(i)*fac_shield(j)*sss enddo C Diagnostics. Comment out or remove after debugging! @@ -5672,7 +5727,7 @@ C This subroutine calculates the excluded-volume interaction energy between C peptide-group centers and side chains and its gradient in virtual-bond and C side-chain vectors. C - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.GEO' include 'COMMON.VAR' @@ -5685,7 +5740,14 @@ C include 'COMMON.CONTROL' include 'COMMON.SPLITELE' integer xshift,yshift,zshift - dimension ggg(3) + double precision ggg(3) + integer i,iint,j,k,iteli,itypj,subchap + double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1, + & fac,e1,e2,rij + double precision evdw2,evdw2_14,evdwij + double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp, + & dist_temp, dist_init + double precision sscale,sscagrad evdw2=0.0D0 evdw2_14=0.0d0 c print *,boxxsize,boxysize,boxzsize,'wymiary pudla' @@ -5831,8 +5893,9 @@ c if (sss.eq.0) print *,'czasem jest OK' endif evdwij=e1+e2 evdw2=evdw2+evdwij*sss - if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') - & 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli), + if (energy_dec) write (iout,'(a6,2i5,3f7.3,2i3,3e11.3)') + & 'evdw2',i,j,1.0d0/dsqrt(rrij),sss, + & evdwij,iteli,itypj,fac,aad(itypj,iteli), & bad(itypj,iteli) C C Calculate contributions to the gradient in the virtual-bond and SC vectors. diff --git a/source/unres/src-HCD-5D/initialize_p.F b/source/unres/src-HCD-5D/initialize_p.F index c73426c..6806e62 100644 --- a/source/unres/src-HCD-5D/initialize_p.F +++ b/source/unres/src-HCD-5D/initialize_p.F @@ -46,6 +46,7 @@ cMS$ATTRIBUTES C :: proc_proc include 'COMMON.DERIV' include 'COMMON.SPLITELE' include 'COMMON.VAR' + include 'COMMON.MD' c Common blocks from the diagonalization routines integer IR,IW,IP,IJK,IPK,IDAF,NAV,IODA,KDIAG,ICORFL,IXDR integer i,idumm,j,k,l,ichir1,ichir2,iblock,m @@ -68,7 +69,7 @@ c NaNQ initialization call proc_proc(rr,i) #endif #endif - + itime_mat=0. kdiag=0 icorfl=0 iw=2 diff --git a/source/unres/src-HCD-5D/lagrangian_lesyng.F b/source/unres/src-HCD-5D/lagrangian_lesyng.F index 4230e10..1180645 100644 --- a/source/unres/src-HCD-5D/lagrangian_lesyng.F +++ b/source/unres/src-HCD-5D/lagrangian_lesyng.F @@ -891,6 +891,7 @@ c--------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.LAGRANGE.5diag' include 'COMMON.INTERACT' + include 'COMMON.VAR' integer ndim double precision forces(3*ndim),accel(3,0:maxres2),rs(ndim), & xsolv(ndim),d_a_vec(6*nres) @@ -904,9 +905,18 @@ Compute accelerations in Calpha and SC innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=iposc,iposc+n-1 - rs(i)=forces(3*(i-1)+j) + rs(i-iposc+1)=forces(3*(i-1)+j) enddo +#ifdef DEBUG + write (iout,*) "j",j," chain",ichain + write (iout,*) "rs" + write (iout,'(f10.5)') (rs(i),i=1,n) +#endif call FDISYS (n,DM(iposc),DU1(iposc),DU2(iposc),rs,xsolv) +#ifdef DEBUG + write (iout,*) "xsolv" + write (iout,'(f10.5)') (xsolv(i),i=1,n) +#endif ind=1 do i=innt,inct if (itype(i).eq.10)then @@ -988,7 +998,7 @@ C Convert d_a to virtual-bon-vector basis enddo #ifdef DEBUG write (iout,*) "d_a_vec" - write (iout,'(3f10.5)') (d_a_vec(j),j=1,dimen3) + write (iout,'(3f10.5)') (d_a_vec(j),j=1,3*(nct-nnt+nside)) #endif return end diff --git a/source/unres/src-HCD-5D/make_xx_list.F b/source/unres/src-HCD-5D/make_xx_list.F new file mode 100644 index 0000000..fb6c055 --- /dev/null +++ b/source/unres/src-HCD-5D/make_xx_list.F @@ -0,0 +1,489 @@ + subroutine make_SCSC_inter_list + implicit none + include "DIMENSIONS" +#ifdef MPI + include 'mpif.h' + include "COMMON.SETUP" +#endif + 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(200*maxres),contlistj(200*maxres) +! integer :: newcontlisti(200*nres),newcontlistj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint, + & ilist_sc,g_ilist_sc + integer displ(0:max_fg_procs),i_ilist_sc(0:max_fg_procs),ierr +! print *,"START make_SC" +#ifdef DEBUG + write (iout,*) "make_SCSC_inter_list" +#endif + r_buff_list=5.0d0 + ilist_sc=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) + xi=dmod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=dmod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=dmod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + 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) + xj=dmod(xj,boxxsize) + if (xj.lt.0) xj=xj+boxxsize + yj=dmod(yj,boxysize) + if (yj.lt.0) yj=yj+boxysize + zj=dmod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize + dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + 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-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +! r_buff_list is a read value for a buffer + if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then +! Here the list is created + ilist_sc=ilist_sc+1 +! this can be substituted by cantor and anti-cantor + contlisti(ilist_sc)=i + contlistj(ilist_sc)=j + + endif + enddo + enddo + enddo +! call MPI_Reduce(ilist_sc,g_ilist_sc,1,& +! MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +! call MPI_Gather(newnss,1,MPI_INTEGER,& +! i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR) +#ifdef MPI +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_sc + do i=1,ilist_sc + write (iout,*) i,contlisti(i),contlistj(i) + enddo +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_sc,g_ilist_sc,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +c write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_sc,1,MPI_INTEGER, + & i_ilist_sc,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_sc(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlisti,ilist_sc,MPI_INTEGER, + & newcontlisti,i_ilist_sc,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Gatherv(contlistj,ilist_sc,MPI_INTEGER, + & newcontlistj,i_ilist_sc,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_sc,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(newcontlisti,g_ilist_sc,MPI_INT,king,FG_COMM, + & IERR) + call MPI_Bcast(newcontlistj,g_ilist_sc,MPI_INT,king,FG_COMM, + & IERR) +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + + else +#endif + g_ilist_sc=ilist_sc + + do i=1,ilist_sc + newcontlisti(i)=contlisti(i) + newcontlistj(i)=contlistj(i) + enddo +#ifdef MPI + endif +#endif +#ifdef DEBUG + write (iout,*) "after GATHERV",g_ilist_sc + do i=1,g_ilist_sc + write (iout,*) i,newcontlisti(i),newcontlistj(i) + enddo +#endif + call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end) + return + end +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine make_SCp_inter_list + implicit none + include "DIMENSIONS" +#ifdef MPI + include 'mpif.h' + include "COMMON.SETUP" +#endif + 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(200*maxres),contlistscpj(200*maxres) +! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint, + & ilist_scp,g_ilist_scp + integer displ(0:max_fg_procs),i_ilist_scp(0:max_fg_procs),ierr + integer contlistscpi_f(200*maxres),contlistscpj_f(200*maxres) + integer ilist_scp_first,ifirstrun,g_ilist_sc +! print *,"START make_SC" +#ifdef DEBUG + write (iout,*) "make_SCp_inter_list" +#endif + r_buff_list=5.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)) + yi=0.5D0*(c(2,i)+c(2,i+1)) + zi=0.5D0*(c(3,i)+c(3,i+1)) + xi=mod(xi,boxxsize) + if (xi.lt.0) xi=xi+boxxsize + yi=mod(yi,boxysize) + if (yi.lt.0) yi=yi+boxysize + zi=mod(zi,boxzsize) + if (zi.lt.0) zi=zi+boxzsize + + 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) + 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-xi)**2+(yj-yi)**2+(zj-zi)**2 + xj_safe=xj + yj_safe=yj + zj_safe=zj + subchap=0 + 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-xi)**2+(yj-yi)**2+(zj-zi)**2 + if(dist_temp.lt.dist_init) then + dist_init=dist_temp + xj_temp=xj + yj_temp=yj + zj_temp=zj + subchap=1 + endif + enddo + enddo + enddo + if (subchap.eq.1) then + xj=xj_temp-xi + yj=yj_temp-yi + zj=zj_temp-zi + else + xj=xj_safe-xi + yj=yj_safe-yi + zj=zj_safe-zi + endif +#ifdef DEBUG + ! r_buff_list is a read value for a buffer + if ((sqrt(dist_init).le.(r_cut_int)).and.(ifirstrun.eq.0)) + & then +! Here the list is created + ilist_scp_first=ilist_scp_first+1 +! this can be substituted by cantor and anti-cantor + contlistscpi_f(ilist_scp_first)=i + contlistscpj_f(ilist_scp_first)=j + endif +#endif +! r_buff_list is a read value for a buffer + if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then +! Here the list is created + ilist_scp=ilist_scp+1 +! this can be substituted by cantor and anti-cantor + contlistscpi(ilist_scp)=i + contlistscpj(ilist_scp)=j + endif + enddo + enddo + enddo +#ifdef MPI +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_scp + do i=1,ilist_scp + write (iout,*) i,contlistscpi(i),contlistscpj(i) + enddo +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_scp,g_ilist_scp,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +c write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_scp,1,MPI_INTEGER, + & i_ilist_scp,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_scp(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistscpi,ilist_scp,MPI_INTEGER, + & newcontlistscpi,i_ilist_scp,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Gatherv(contlistscpj,ilist_scp,MPI_INTEGER, + & newcontlistscpj,i_ilist_scp,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_scp,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(newcontlistscpi,g_ilist_scp,MPI_INT,king,FG_COMM, + & IERR) + call MPI_Bcast(newcontlistscpj,g_ilist_scp,MPI_INT,king,FG_COMM, + & IERR) +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + else +#endif + g_ilist_scp=ilist_scp + + do i=1,ilist_scp + newcontlistscpi(i)=contlistscpi(i) + newcontlistscpj(i)=contlistscpj(i) + enddo +#ifdef MPI + endif +#endif +#ifdef DEBUG + write (iout,*) "after MPIREDUCE",g_ilist_scp + do i=1,g_ilist_scp + write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i) + enddo + +! if (ifirstrun.eq.0) ifirstrun=1 +! do i=1,ilist_scp_first +! do j=1,g_ilist_scp +! if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.& +! (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126 +! enddo +! print *,itime_mat,"ERROR matrix needs updating" +! print *,contlistscpi_f(i),contlistscpj_f(i) +! 126 continue +! enddo +#endif + call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end) + + return + end +!----------------------------------------------------------------------------- + subroutine make_pp_inter_list + implicit none + include "DIMENSIONS" +#ifdef MPI + include 'mpif.h' + include "COMMON.SETUP" +#endif + 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 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 contlistppi(200*maxres),contlistppj(200*maxres) +! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres) + integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint, + & ilist_pp,g_ilist_pp + integer displ(0:max_fg_procs),i_ilist_pp(0:max_fg_procs),ierr +! print *,"START make_SC" +#ifdef DEBUG + write (iout,*) "make_pp_inter_list" +#endif + ilist_pp=0 + r_buff_list=5.0 + do i=iatel_s,iatel_e + 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 + do j=ielstart(i),ielend(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 + + if (sqrt(dist_init).le.(r_cut_int+r_buff_list)) then +! Here the list is created + ilist_pp=ilist_pp+1 +! this can be substituted by cantor and anti-cantor + contlistppi(ilist_pp)=i + contlistppj(ilist_pp)=j + endif + enddo + enddo +! enddo +#ifdef MPI +#ifdef DEBUG + write (iout,*) "before MPIREDUCE",ilist_pp + do i=1,ilist_pp + write (iout,*) i,contlistppi(i),contlistppj(i) + enddo +#endif + if (nfgtasks.gt.1)then + + call MPI_Reduce(ilist_pp,g_ilist_pp,1, + & MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR) +! write(iout,*) "before bcast",g_ilist_sc + call MPI_Gather(ilist_pp,1,MPI_INTEGER, + & i_ilist_pp,1,MPI_INTEGER,king,FG_COMM,IERR) + displ(0)=0 + do i=1,nfgtasks-1,1 + displ(i)=i_ilist_pp(i-1)+displ(i-1) + enddo +! write(iout,*) "before gather",displ(0),displ(1) + call MPI_Gatherv(contlistppi,ilist_pp,MPI_INTEGER, + & newcontlistppi,i_ilist_pp,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Gatherv(contlistppj,ilist_pp,MPI_INTEGER, + & newcontlistppj,i_ilist_pp,displ,MPI_INTEGER, + & king,FG_COMM,IERR) + call MPI_Bcast(g_ilist_pp,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(newcontlistppi,g_ilist_pp,MPI_INT,king,FG_COMM, + & IERR) + call MPI_Bcast(newcontlistppj,g_ilist_pp,MPI_INT,king,FG_COMM, + & IERR) + +! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM) + + else +#endif + g_ilist_pp=ilist_pp + + do i=1,ilist_pp + newcontlistppi(i)=contlistppi(i) + newcontlistppj(i)=contlistppj(i) + enddo +#ifdef MPI + endif +#endif + call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end) +#ifdef DEBUG + write (iout,*) "after MPIREDUCE",g_ilist_pp + do i=1,g_ilist_pp + write (iout,*) i,newcontlistppi(i),newcontlistppj(i) + enddo +#endif + return + end diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 4765a41..8f4d05d 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -171,6 +171,7 @@ c call readi(controlcard,'IZ_SC',iz_sc,0) timlim=60.0D0*timlim safety = 60.0d0*safety modecalc=0 + call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100) call reada(controlcard,"T_BATH",t_bath,300.0d0) minim=(index(controlcard,'MINIMIZE').gt.0) dccart=(index(controlcard,'CART').gt.0) @@ -546,6 +547,10 @@ c if performing umbrella sampling, fragments constrained are read from the frag & "Initial time step of numerical integration:",d_time, & " natural units" write (iout,'(60x,f10.5,a)') d_time*48.9," fs" + write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int, + & " A" + write(iout,'(a60,i5)')"Frequency of updating interaction list", + & imatupdate if (RESPA) then write (iout,'(2a,i4,a)') & "A-MTS algorithm used; initial time step for fast-varying", diff --git a/source/unres/src-HCD-5D/stochfric.F b/source/unres/src-HCD-5D/stochfric.F index b8069d9..09b6877 100644 --- a/source/unres/src-HCD-5D/stochfric.F +++ b/source/unres/src-HCD-5D/stochfric.F @@ -70,6 +70,9 @@ c write (iout,*) "friction_force j",j," ichain",ichain, c & " n",n," iposc",iposc,iposc+n-1 innt=chain_border(1,ichain) inct=chain_border(2,ichain) +c diagnostics +c innt=chain_border(1,1) +c inct=chain_border(2,1) do i=innt,inct vvec(ind+1)=v_work(j,i) ind=ind+1 @@ -324,6 +327,10 @@ c Compute the stochastic forces acting on bodies. Store in force. innt=chain_border(1,ichain) inct=chain_border(2,ichain) iposc=iposd_chain(ichain) +c for debugging only +c innt=chain_border(1,1) +c inct=chain_border(2,1) +c iposc=iposd_chain(1) c write (iout,*)"stochastic_force ichain=",ichain," innt",innt, c & " inct",inct," iposc",iposc do j=1,3 -- 1.7.9.5