& 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
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)
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
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)
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
#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
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
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)
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
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.
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
& sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision boxshift
+ double precision x0,y0,r012,rij12,facx0,
+ & facx02,afacx0,bfacx0,abfacx0,Afac,BBfac,Afacsig,Bfacsig
+c alpha_GB=0.5d0
+c alpha_GB=0.01d0
+c alpha_GB1=1.0d0+1.0d0/alpha_GB
evdw=0.0D0
ccccc energy_dec=.false.
C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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)')
& 'evdw',i,j,evdwij,' ss'
C triple bond artifac removal
- do k=j+1,iend(i,iint)
+c do k=j+1,iend(i,iint)
+ do k=j+1,nct
C search over all next residues
if (dyn_ss_mask(k)) then
C check if they are cysteins
c for diagnostics; uncomment
c rij_shift=1.2*sig0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
- if (rij_shift.le.0.0D0) then
- evdw=1.0D20
+c if (rij_shift.le.0.0D0) then
+ x0=alpha_GB*(sig-sig0ij)
+ if (energy_dec) write (iout,*) i,j," x0",x0
+ if (rij_shift.le.x0) then
+c sig=2.0d0*sig0ij
+ sigder=-sig*sigsq
+c sigder=0.0d0
+ fac=rij**expon
+ rij12=fac*fac
+c rij12=1.0d0
+ x0=alpha_GB*(sig-sig0ij)
+ facx0=1.0d0/x0**expon
+ facx02=facx0*facx0
+ r012=((1.0d0+alpha_GB)*(sig-sig0ij))**(2*expon)
+ afacx0=aa*facx02
+ bfacx0=bb*facx0
+ abfacx0=afacx0+0.5d0*bfacx0
+ Afac=alpha_GB1*abfacx0
+ Afacsig=0.5d0*alpha_GB1*bfacx0/(sig-sig0ij)
+ BBfac=Afac-(afacx0+bfacx0)
+c BBfac=0.0d0
+ Bfacsig=(-alpha_GB1*(abfacx0+afacx0)+
+ & (afacx0+afacx0+bfacx0))/(sig-sig0ij)
+c Bfacsig=0.0d0
+ Afac=Afac*r012
+ Afacsig=Afacsig*r012
+c Afac=1.0d0
+c Afacsig=0.0d0
+c w(x)=4*eps*((1.0+1.0/alpha_GB)*(y0**12-0.5*y0**6)*(r0/x)**12-(1+1/alpha)*(y0**12-0.5*y0**6)+y0**12-y0**6)
+c eps1=1.0d0
+c eps2rt=1.0d0
+c eps3rt=1.0d0
+ e1 = eps1*eps2rt*eps3rt*Afac*rij12
+ e2 = -eps1*eps2rt*eps3rt*BBfac
+ evdwij = e1+e2
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+c eps2der=0.0d0
+c eps3der=0.0d0
+c eps1_om12=0.0d0
+ evdwij=evdwij*eps2rt*eps3rt
+c Afacsig=0.d0
+c Bfacsig=0.0d0
+ if (lprn) then
+ write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),18(0pf9.5))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & eps1*eps2rt**2*eps3rt**2,om1,om2,om12,
+ & 1.0D0/rij,rij_shift,
+ & evdwij
+ endif
+ if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
+ & 'RE r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+ evdw=evdw+evdwij*sss
+C Calculate gradient components.
+ e1=e1*eps2rt*eps3rt
+ sigder=-expon*eps1*eps2rt*eps2rt*eps3rt*eps3rt
+ & *(Afacsig*rij12-Bfacsig)*sigder
+ fac=-2.0d0*expon*e1*rij*rij
+c print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c & evdwij,fac,sigma(itypi,itypj),expon
+ fac=fac+evdwij*sssgrad/sss*rij
+c fac=0.0d0
+c write (iout,*) "sigder",sigder," fac",fac," e1",e1,
+c & " e2",e2," sss",sss," sssgrad",sssgrad,"esp123",
+c & eps1*eps2rt**2*eps3rt**2
+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)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+C gg_lipi(3)=0.0d0
+C gg_lipj(3)=0.0d0
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+c evdw=1.0D20
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
- endif
- sigder=-sig*sigsq
+ else
+ rij_shift=1.0D0/rij_shift
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
+ fac=rij_shift**expon
C here to start with
C if (c(i,3).gt.
- faclip=fac
- e1=fac*fac*aa
- e2=fac*bb
- evdwij=eps1*eps2rt*eps3rt*(e1+e2)
- eps2der=evdwij*eps3rt
- eps3der=evdwij*eps2rt
+ faclip=fac
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
C &((sslipi+sslipj)/2.0d0+
C &(2.0d0-sslipi-sslipj)/2.0d0)
c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
- evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
- if (lprn) then
- sigm=dabs(aa/bb)**(1.0D0/6.0D0)
- epsi=bb**2/aa
- write (iout,'(2(a3,i3,2x),17(0pf7.3))')
- & restyp(itypi),i,restyp(itypj),j,
- & epsi,sigm,chi1,chi2,chip1,chip2,
- & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
- & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
- & evdwij
- endif
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
+ & 'GB r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
+ if (lprn) then
+ write (iout,*) "aa",aa," bb",bb," sig0ij",sig0ij
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij
+ endif
- if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
- & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
C Calculate gradient components.
- e1=e1*eps1*eps2rt**2*eps3rt**2
- fac=-expon*(e1+evdwij)*rij_shift
- sigder=fac*sigder
- fac=rij*fac
-c print '(2i4,6f8.4)',i,j,sss,sssgrad*
-c & evdwij,fac,sigma(itypi,itypj),expon
- fac=fac+evdwij*sssgrad/sss*rij
-c fac=0.0d0
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
+c print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c & evdwij,fac,sigma(itypi,itypj),expon
+ fac=fac+evdwij*sssgrad/sss*rij
+c fac=0.0d0
C Calculate the radial part of the gradient
- gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ 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)))
- gg_lipj(3)=ssgradlipj*gg_lipi(3)
- gg_lipi(3)=gg_lipi(3)*ssgradlipi
-C gg_lipi(3)=0.0d0
-C gg_lipj(3)=0.0d0
- gg(1)=xj*fac
- gg(2)=yj*fac
- gg(3)=zj*fac
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+C gg_lipi(3)=0.0d0
+C gg_lipj(3)=0.0d0
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
+ endif
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
call sc_grad
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
C side-chain vectors.
C
implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
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,
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'
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
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)
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)
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
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
projects / unres.git / blobdiff