!
implicit none
!-----------------------------------------------------------------------------
+! Max. number of contacts per residue
+! integer :: maxconts
+!-----------------------------------------------------------------------------
+! Max. number of derivatives of virtual-bond and side-chain vectors in theta
+! or phi.
+! integer :: maxdim
+!-----------------------------------------------------------------------------
+! Max. number of SC contacts
+! integer :: maxcont
+!-----------------------------------------------------------------------------
! Max. number of variables
integer :: maxvar
!-----------------------------------------------------------------------------
-! Max number of torsional terms in SCCOR
- integer,parameter :: maxterm_sccor=6
+! Max number of torsional terms in SCCOR in control_data
+! integer,parameter :: maxterm_sccor=6
!-----------------------------------------------------------------------------
! Maximum number of SC local term fitting function coefficiants
integer,parameter :: maxsccoef=65
!-----------------------------------------------------------------------------
+! commom.calc common/calc/
+!-----------------------------------------------------------------------------
! commom.contacts
! common /contacts/
! Change 12/1/95 - common block CONTACTS1 included.
subroutine etotal(energia)
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
- use MD_data, only: totT
+ use MD_data
#ifndef ISNAN
external proc_proc
#ifdef WINPGI
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
+! shielding effect varibles for MPI
+! real(kind=8) fac_shieldbuf(maxres),
+! & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
+! & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
+! & grad_shieldbuf(3,-1:maxres)
+! integer ishield_listbuf(maxres),
+! &shield_listbuf(maxcontsshi,maxres)
+
! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
! & " nfgtasks",nfgtasks
if (nfgtasks.gt.1) then
!
! Compute the side-chain and electrostatic interaction energy
!
- goto (101,102,103,104,105,106) ipot
+! goto (101,102,103,104,105,106) ipot
+ select case(ipot)
! Lennard-Jones potential.
- 101 call elj(evdw)
+! 101 call elj(evdw)
+ case (1)
+ call elj(evdw)
!d print '(a)','Exit ELJcall el'
- goto 107
+! goto 107
! Lennard-Jones-Kihara potential (shifted).
- 102 call eljk(evdw)
- goto 107
+! 102 call eljk(evdw)
+ case (2)
+ call eljk(evdw)
+! goto 107
! Berne-Pechukas potential (dilated LJ, angular dependence).
- 103 call ebp(evdw)
- goto 107
+! 103 call ebp(evdw)
+ case (3)
+ call ebp(evdw)
+! goto 107
! Gay-Berne potential (shifted LJ, angular dependence).
- 104 call egb(evdw)
- goto 107
+! 104 call egb(evdw)
+ case (4)
+ call egb(evdw)
+! goto 107
! Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
- 105 call egbv(evdw)
- goto 107
+! 105 call egbv(evdw)
+ case (5)
+ call egbv(evdw)
+! goto 107
! Soft-sphere potential
- 106 call e_softsphere(evdw)
+! 106 call e_softsphere(evdw)
+ case (6)
+ call e_softsphere(evdw)
!
! Calculate electrostatic (H-bonding) energy of the main chain.
!
- 107 continue
+! 107 continue
+ case default
+ write(iout,*)"Wrong ipot"
+! return
+! 50 continue
+ end select
+! continue
!mc
!mc Sep-06: egb takes care of dynamic ss bonds too
.or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
#endif
call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+! write (iout,*) "ELEC calc"
else
ees=0.0d0
evdw1=0.0d0
! Calculate excluded-volume interaction energy between peptide groups
! and side chains.
!
+!elwrite(iout,*) "in etotal calc exc;luded",ipot
if (ipot.lt.6) then
if(wscp.gt.0d0) then
! write (iout,*) "Soft-sphere SCP potential"
call escp_soft_sphere(evdw2,evdw2_14)
endif
+!elwrite(iout,*) "in etotal before ebond",ipot
!
! Calculate the bond-stretching energy
!
call ebond(estr)
+!elwrite(iout,*) "in etotal afer ebond",ipot
+
!
! Calculate the disulfide-bridge and other energy and the contributions
! from other distance constraints.
- print *,'Calling EHPB'
+! print *,'Calling EHPB'
call edis(ehpb)
+!elwrite(iout,*) "in etotal afer edis",ipot
! print *,'EHPB exitted succesfully.'
!
! Calculate the virtual-bond-angle energy.
! Calculate the SC local energy.
!
call esc(escloc)
+!elwrite(iout,*) "in etotal afer esc",ipot
! print *,"Processor",myrank," computed USC"
!
! Calculate the virtual-bond torsional energy.
!
! 6/23/01 Calculate double-torsional energy
!
+!elwrite(iout,*) "in etotal",ipot
if (wtor_d.gt.0) then
call etor_d(etors_d)
else
ecorr6=0.0d0
eturn6=0.0d0
endif
+!elwrite(iout,*) "in etotal",ipot
if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then
call multibody_hb(ecorr,ecorr5,ecorr6,n_corr,n_corr1)
!d write (iout,*) "multibody_hb ecorr",ecorr
endif
+!elwrite(iout,*) "afeter multibody hb"
! print *,"Processor",myrank," computed Ucorr"
!
! If performing constraint dynamics, call the constraint energy
! after the equilibration time
if(usampl.and.totT.gt.eq_time) then
+!elwrite(iout,*) "afeter multibody hb"
call EconstrQ
+!elwrite(iout,*) "afeter multibody hb"
call Econstr_back
+!elwrite(iout,*) "afeter multibody hb"
else
Uconst=0.0d0
Uconst_back=0.0d0
endif
+!elwrite(iout,*) "after Econstr"
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#ifdef TIMING
time_sumene=time_sumene+MPI_Wtime()-time00
#endif
+!el call enerprint(energia)
+!elwrite(iout,*)"finish etotal"
return
end subroutine etotal
!-----------------------------------------------------------------------------
integer :: ierr
real(kind=8) :: time00
if (nfgtasks.gt.1 .and. reduce) then
-!el #define DEBUG
+
#ifdef DEBUG
write (iout,*) "energies before REDUCE"
call enerprint(energia)
#ifdef MPI
endif
#endif
-!el #undef DUBUG
+! call enerprint(energia)
call flush(iout)
return
end subroutine sum_energy
real(kind=8) :: kfac=2.4d0
real(kind=8) :: x,x2,x3,x4,x5,licznik=1.12692801104297249644
!el local variables
- real(kind=8) :: t_bath,facT,facT2,facT3,facT4,facT5
+ real(kind=8) :: t_bath,facT(6) !,facT2,facT3,facT4,facT5,facT6
+ real(kind=8) :: T0=3.0d2
integer :: ierror
! facT=temp0/t_bath
! facT=2*temp0/(t_bath+temp0)
if (rescale_mode.eq.0) then
- facT=1.0d0
- facT2=1.0d0
- facT3=1.0d0
- facT4=1.0d0
- facT5=1.0d0
+ facT(1)=1.0d0
+ facT(2)=1.0d0
+ facT(3)=1.0d0
+ facT(4)=1.0d0
+ facT(5)=1.0d0
+ facT(6)=1.0d0
else if (rescale_mode.eq.1) then
- facT=kfac/(kfac-1.0d0+t_bath/temp0)
- facT2=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2)
- facT3=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3)
- facT4=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4)
- facT5=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5)
+ facT(1)=kfac/(kfac-1.0d0+t_bath/temp0)
+ facT(2)=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2)
+ facT(3)=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3)
+ facT(4)=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4)
+ facT(5)=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5)
+#ifdef WHAM_RUN
+!#if defined(WHAM_RUN) || defined(CLUSTER)
+#if defined(FUNCTH)
+! tt = 1.0d0/(beta_h(ib,ipar)*1.987D-3)
+ facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0
+#elif defined(FUNCT)
+ facT(6)=t_bath/T0
+#else
+ facT(6)=1.0d0
+#endif
+#endif
else if (rescale_mode.eq.2) then
x=t_bath/temp0
x2=x*x
x3=x2*x
x4=x3*x
x5=x4*x
- facT=licznik/dlog(dexp(x)+dexp(-x))
- facT2=licznik/dlog(dexp(x2)+dexp(-x2))
- facT3=licznik/dlog(dexp(x3)+dexp(-x3))
- facT4=licznik/dlog(dexp(x4)+dexp(-x4))
- facT5=licznik/dlog(dexp(x5)+dexp(-x5))
+ facT(1)=licznik/dlog(dexp(x)+dexp(-x))
+ facT(2)=licznik/dlog(dexp(x2)+dexp(-x2))
+ facT(3)=licznik/dlog(dexp(x3)+dexp(-x3))
+ facT(4)=licznik/dlog(dexp(x4)+dexp(-x4))
+ facT(5)=licznik/dlog(dexp(x5)+dexp(-x5))
+#ifdef WHAM_RUN
+!#if defined(WHAM_RUN) || defined(CLUSTER)
+#if defined(FUNCTH)
+ facT(6)=(320.0+80.0*dtanh((t_bath-320.0)/80.0))/320.0
+#elif defined(FUNCT)
+ facT(6)=t_bath/T0
+#else
+ facT(6)=1.0d0
+#endif
+#endif
else
write (iout,*) "Wrong RESCALE_MODE",rescale_mode
write (*,*) "Wrong RESCALE_MODE",rescale_mode
#endif
stop 555
endif
- welec=weights(3)*fact
- wcorr=weights(4)*fact3
- wcorr5=weights(5)*fact4
- wcorr6=weights(6)*fact5
- wel_loc=weights(7)*fact2
- wturn3=weights(8)*fact2
- wturn4=weights(9)*fact3
- wturn6=weights(10)*fact5
- wtor=weights(13)*fact
- wtor_d=weights(14)*fact2
- wsccor=weights(21)*fact
+ welec=weights(3)*fact(1)
+ wcorr=weights(4)*fact(3)
+ wcorr5=weights(5)*fact(4)
+ wcorr6=weights(6)*fact(5)
+ wel_loc=weights(7)*fact(2)
+ wturn3=weights(8)*fact(2)
+ wturn4=weights(9)*fact(3)
+ wturn6=weights(10)*fact(5)
+ wtor=weights(13)*fact(1)
+ wtor_d=weights(14)*fact(2)
+ wsccor=weights(21)*fact(1)
return
end subroutine rescale_weights
! include 'COMMON.SBRIDGE'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
-
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer :: ii
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
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
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
evdw=evdw+evdwij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') &
'evdw',i,j,evdwij,' ss'
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,' ss'
ELSE
!el ind=ind+1
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
! dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
-! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
-! & 1.0d0/vbld(j+nres)
+! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,&
+! 1.0d0/vbld(j+nres) !d
! write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
! alf1=0.0D0
! alf2=0.0D0
! alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ 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
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
! write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
-! write (iout,*) "j",j," dc_norm",
-! & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+! write (iout,*) "j",j," dc_norm",& !d
+! dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+! write(iout,*)"rrij ",rrij
+! write(iout,*)"xj yj zj ", xj, yj, zj
+! write(iout,*)"xi yi zi ", xi, yi, zi
+! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
! Calculate angle-dependent terms of energy and contributions to their
! derivatives.
call sc_angular
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
rij_shift=1.0D0/rij-sig+sig0ij
+! write(iout,*)" rij_shift",rij_shift," rij",rij," sig",sig,&
+! "sig0ij",sig0ij
! for diagnostics; uncomment
! rij_shift=1.2*sig0ij
! I hate to put IF's in the loops, but here don't have another choice!!!!
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
-! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,&
-! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+! write(iout,*)"aa, bb ",aa(:,:),bb(:,:)
+! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
+! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
endif
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
- 'evdw',i,j,evdwij
+ 'evdw',i,j,evdwij !,"egb"
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+! print *,'before fac',fac,rij,evdwij
+ fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij
+! print *,'grad part scale',fac, &
+! evdwij*sss_ele_grad/sss_ele_cut &
+! /sigma(itypi,itypj)*rij
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+! print *,'before sc_grad', gg(1),gg(2),gg(3)
! Calculate angular part of the gradient.
call sc_grad
ENDIF ! dyn_ss
integer :: i,iti1,iti,k,l
real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2
-! allocate(Ug(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2(2,2,nres)) !(2,2,maxres)
-! allocate(Ugder(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2der(2,2,nres)) !(2,2,maxres)
-! allocate(obrot(2,nres)) !(2,maxres)
-! allocate(obrot2(2,nres)) !(2,maxres)
-! allocate(obrot_der(2,nres)) !(2,maxres)
-! allocate(obrot2_der(2,nres)) !(2,maxres)
-! allocate(costab2(nres)) !(maxres)
-! allocate(sintab2(nres)) !(maxres)
-! allocate(costab(nres)) !(maxres)
-! allocate(sintab(nres)) !(maxres)
-
-! allocate(Ub2(2,nres)) !(2,maxres)
-! allocate(Ctobr(2,nres)) !(2,maxres)
-! allocate(Dtobr2(2,nres)) !(2,maxres)
-! allocate(mu(2,nres)) !(2,maxres)
-! allocate(muder(2,nres)) !(2,maxres)
-! allocate(Ub2der(2,nres)) !(2,maxres)
-! allocate(Ctobrder(2,nres)) !(2,maxres)
-! allocate(Dtobr2der(2,nres)) !(2,maxres)
-
-! allocate(EUg(2,2,nres)) !(2,2,maxres)
-! allocate(CUg(2,2,nres)) !(2,2,maxres)
-! allocate(DUg(2,2,nres)) !(2,2,maxres)
-! allocate(DtUg2(2,2,nres)) !(2,2,maxres)
-! allocate(EUgder(2,2,nres)) !(2,2,maxres)
-! allocate(CUgder(2,2,nres)) !(2,2,maxres)
-! allocate(DUgder(2,2,nres)) !(2,2,maxres)
-! allocate(Dtug2der(2,2,nres)) !(2,2,maxres)
-
-! allocate(Ug2Db1t(2,nres)) !(2,maxres)
-! allocate(Ug2Db1tder(2,nres)) !(2,maxres)
-! allocate(CUgb2(2,nres)) !(2,maxres)
-! allocate(CUgb2der(2,nres)) !(2,maxres)
-
-! allocate(EUgC(2,2,nres)) !(2,2,maxres)
-! allocate(EUgCder(2,2,nres)) !(2,2,maxres)
-! allocate(EUgD(2,2,nres)) !(2,2,maxres)
-! allocate(EUgDder(2,2,nres)) !(2,2,maxres)
-! allocate(DtUg2EUg(2,2,nres)) !(2,2,maxres)
-! allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres)
-
-! allocate(Ug2DtEUgder(2,2,2,nres)) !(2,2,2,maxres)
-! allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres)
-
!
! Compute the virtual-bond-torsional-angle dependent quantities needed
! to calculate the el-loc multibody terms of various order.
!
+!AL el mu=0.0d0
#ifdef PARMAT
do i=ivec_start+2,ivec_end+2
#else
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
enddo
-!d write (iout,*) 'mu ',mu(:,i-2)
+! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2)
+! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1)
+! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2)
!d write (iout,*) 'mu1',mu1(:,i-2)
!d write (iout,*) 'mu2',mu2(:,i-2)
if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
enddo
endif
#if defined(MPI) && defined(PARMAT)
-!el #define DUBUG
#ifdef DEBUG
! if (fg_rank.eq.0) then
write (iout,*) "Arrays UG and UGDER before GATHER"
!d enddo
!d enddo
return
-!el #undef DUBUG
end subroutine set_matrices
!-----------------------------------------------------------------------------
subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
!d enddo
!d call check_vecgrad
!d stop
+! ees=0.0d0 !AS
+! evdw1=0.0d0
+! eel_loc=0.0d0
+! eello_turn3=0.0d0
+! eello_turn4=0.0d0
+ t_eelecij=0.0d0
+ ees=0.0D0
+ evdw1=0.0D0
+ eel_loc=0.0d0
+ eello_turn3=0.0d0
+ eello_turn4=0.0d0
+!
+
if (icheckgrad.eq.1) then
+!el
+! do i=0,2*nres+2
+! dc_norm(1,i)=0.0d0
+! dc_norm(2,i)=0.0d0
+! dc_norm(3,i)=0.0d0
+! enddo
do i=1,nres-1
fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
do k=1,3
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
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
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
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
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
+
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
! include 'COMMON.VECTORS'
! include 'COMMON.FFIELD'
! include 'COMMON.TIME1'
- real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg
+ real(kind=8),dimension(3) :: ggg,gggp,gggm,erij,dcosb,dcosg,xtemp
real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
!el integer :: num_conti,j1,j2
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el dz_normi,xmedi,ymedi,zmedi
0.0d0,0.0d0,1.0d0/),shape(unmat))
! integer :: maxconts=nres/4
!el local variables
- integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+ integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
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-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
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+!C print *,i,j
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
+
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0D0/rij
rij=dsqrt(rij)
+!C print *,xmedi,ymedi,zmedi,xj,yj,zj,boxxsize,rij
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! sss_ele_cut=1.0d0
+! sss_ele_grad=0.0d0
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+! if (sss_ele_cut.le.0.0) go to 128
+
rmij=1.0D0/rij
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
eesij=el1+el2
! 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
- ees=ees+eesij
- evdw1=evdw1+evdwij
+ ees=ees+eesij*sss_ele_cut
+ evdw1=evdw1+evdwij*sss_ele_cut
!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
!d & xmedi,ymedi,zmedi,xj,yj,zj
if (energy_dec) then
- write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &
- 'evdw1',i,j,evdwij,&
- iteli,itelj,aaa,evdw1
+! write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)') &
+! 'evdw1',i,j,evdwij,&
+! iteli,itelj,aaa,evdw1
+ write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
endif
!
! Calculate contributions to the Cartesian gradient.
!
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
- facel=-3*rrmij*(el1+eesij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
+ ggg(1)=facel*xj+sss_ele_grad*rmij*eesij*xj
+ ggg(2)=facel*yj+sss_ele_grad*rmij*eesij*yj
+ ggg(3)=facel*zj+sss_ele_grad*rmij*eesij*zj
+
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!grad gelc(l,k)=gelc(l,k)+ggg(l)
!grad enddo
!grad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
!grad enddo
!grad enddo
#else
- facvdw=ev1+evdwij
- facel=el1+eesij
+ facvdw=(ev1+evdwij)*sss_ele_cut
+ facel=(el1+eesij)*sss_ele_cut
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+ ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+ ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
!d print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
!d & (dcosg(k),k=1,3)
do k=1,3
- ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss_ele_cut
enddo
! do k=1,3
! ghalf=0.5D0*ggg(k)
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_ele_cut
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_ele_cut
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+
IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 &
.or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 &
.or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
! Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
-!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! eel_loc_ij=0.0
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
+! if (energy_dec) write (iout,*) "a22",a22," a23",a23," a32",a32," a33",a33
+! if (energy_dec) write (iout,*) "muij",muij
! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
-
- eel_loc=eel_loc+eel_loc_ij
+
+ eel_loc=eel_loc+eel_loc_ij*sss_ele_cut
! Partial derivatives in virtual-bond dihedral angles gamma
if (i.gt.1) &
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)
+ (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)) &
+ *sss_ele_cut
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)
+ (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)) &
+ *sss_ele_cut
! Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
- do l=1,3
- ggg(l)=agg(l,1)*muij(1)+ &
- agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)
+! do l=1,3
+! ggg(1)=(agg(1,1)*muij(1)+ &
+! agg(1,2)*muij(2)+agg(1,3)*muij(3)+agg(1,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*xj
+! ggg(2)=(agg(2,1)*muij(1)+ &
+! agg(2,2)*muij(2)+agg(2,3)*muij(3)+agg(2,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*yj
+! ggg(3)=(agg(3,1)*muij(1)+ &
+! agg(3,2)*muij(2)+agg(3,3)*muij(3)+agg(3,4)*muij(4)) &
+! *sss_ele_cut &
+! +eel_loc_ij*sss_ele_grad*rmij*zj
+ xtemp(1)=xj
+ xtemp(2)=yj
+ xtemp(3)=zj
+
+ do l=1,3
+ ggg(l)=(agg(l,1)*muij(1)+ &
+ agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))&
+ *sss_ele_cut &
+ +eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+
gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
!grad ghalf=0.5d0*ggg(l)
!grad enddo
! 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)
- 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)
- 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)
- 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)
+ 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))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
+ 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))&
+ *sss_ele_cut
+!+eel_loc_ij*sss_ele_grad*rmij*xtemp(l)
enddo
ENDIF
! Change 12/26/95 to calculate four-body contributions to H-bonding energy
r0ij=2.20D0*rpp(iteli,itelj)
! r0ij=1.55D0*rpp(iteli,itelj)
call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
+!elwrite(iout,*) "num_conti",num_conti, "maxconts",maxconts
if (fcont.gt.0.0D0) then
num_conti=num_conti+1
if (num_conti.gt.maxconts) then
ees0mij=0
endif
! ees0mij=0.0D0
- ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij) &
+ *sss_ele_cut
+
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij) &
+ *sss_ele_cut
+
! Diagnostics. Comment out or remove after debugging!
! ees0p(num_conti,i)=0.5D0*fac3*ees0pij
! ees0m(num_conti,i)=0.5D0*fac3*ees0mij
gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
enddo
- gggp(1)=gggp(1)+ees0pijp*xj
- gggp(2)=gggp(2)+ees0pijp*yj
- gggp(3)=gggp(3)+ees0pijp*zj
- gggm(1)=gggm(1)+ees0mijp*xj
- gggm(2)=gggm(2)+ees0mijp*yj
- gggm(3)=gggm(3)+ees0mijp*zj
+ gggp(1)=gggp(1)+ees0pijp*xj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+ gggp(2)=gggp(2)+ees0pijp*yj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+ gggp(3)=gggp(3)+ees0pijp*zj &
+ +ees0p(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
+ gggm(1)=gggm(1)+ees0mijp*xj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*xj*sss_ele_grad
+
+ gggm(2)=gggm(2)+ees0mijp*yj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*yj*sss_ele_grad
+
+ gggm(3)=gggm(3)+ees0mijp*zj &
+ +ees0m(num_conti,i)/sss_ele_cut*rmij*zj*sss_ele_grad
+
! Derivatives due to the contact function
gacont_hbr(1,num_conti,i)=fprimcont*xj
gacont_hbr(2,num_conti,i)=fprimcont*yj
!grad 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)
+ + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut
+
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)
- gacontp_hb3(k,num_conti,i)=gggp(k)
+ + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)&
+ *sss_ele_cut
+
+ gacontp_hb3(k,num_conti,i)=gggp(k) &
+ *sss_ele_cut
+
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)
+ + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) &
+ *sss_ele_cut
+
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)
- gacontm_hb3(k,num_conti,i)=gggm(k)
+ + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) &
+ *sss_ele_cut
+
+ gacontm_hb3(k,num_conti,i)=gggm(k) &
+ *sss_ele_cut
+
enddo
! Diagnostics. Comment out or remove after debugging!
!diag do k=1,3
enddo
endif
endif
+ 128 continue
! t_eelecij=t_eelecij+MPI_Wtime()-time00
return
end subroutine eelecij
! include 'COMMON.CONTROL'
real(kind=8),dimension(3) :: ggg
!el local variables
- integer :: i,iint,j,k,iteli,itypj
+ integer :: i,iint,j,k,iteli,itypj,subchap
real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
- e1,e2,evdwij
+ e1,e2,evdwij,rij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
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)
! 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)-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
+
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(1.0d0/rrij)
+ sss_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! (rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+e1+e2
+ evdw2_14=evdw2_14+(e1+e2)*sss_ele_cut
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)') &
- 'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),&
- bad(itypj,iteli)
+ evdw2=evdw2+evdwij*sss_ele_cut
+! 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,0pf7.3)') &
+ 'evdw2',i,j,evdwij
!
! Calculate contributions to the gradient in the virtual-bond and SC vectors.
!
- fac=-(evdwij+e1)*rrij
+ fac=-(evdwij+e1)*rrij*sss_ele_cut
+ fac=fac+evdwij*sss_ele_grad/rij/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
do i=ibondp_start,ibondp_end
+ if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
- *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*) &
- "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C do j=1,3
+!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C *dc(j,i-1)/vbld(i)
+!C enddo
+!C if (energy_dec) write(iout,*) &
+!C "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+ diff = vbld(i)-vbldpDUM
else
diff = vbld(i)-vbldp0
- if (energy_dec) write (iout,*) &
+ endif
+ if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
"estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
do j=1,3
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
! write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
- endif
+! endif
enddo
estr=0.5d0*AKP*estr+estr1
!
if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
#ifdef OSF
- phii=phi(i)
+ phii=phi(i)
if (phii.ne.phii) phii=150.0
#else
phii=phi(i)
endif
if (i.lt.nres .and. itype(i).ne.ntyp1) then
#ifdef OSF
- phii1=phi(i+1)
+ phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
phii1=pinorm(phii1)
z(1)=cos(phii1)
etheta=0.0D0
do i=ithet_start,ithet_end
if (itype(i-1).eq.ntyp1) cycle
+ if (itype(i-2).eq.ntyp1.or.itype(i).eq.ntyp1) cycle
if (iabs(itype(i+1)).eq.20) iblock=2
if (iabs(itype(i+1)).ne.20) iblock=1
dethetai=0.0d0
coskt(k)=dcos(k*theti2)
sinkt(k)=dsin(k*theti2)
enddo
- if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+ if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
#ifdef OSF
phii=phi(i)
if (phii.ne.phii) phii=150.0
enddo
else
phii=0.0d0
- ityp1=nthetyp+1
+ ityp1=ithetyp(itype(i-2))
do k=1,nsingle
cosph1(k)=0.0d0
sinph1(k)=0.0d0
enddo
endif
- if (i.lt.nres .and. itype(i).ne.ntyp1) then
+ if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
#ifdef OSF
phii1=phi(i+1)
if (phii1.ne.phii1) phii1=150.0
enddo
else
phii1=0.0d0
- ityp3=nthetyp+1
+ ityp3=ithetyp(itype(i))
do k=1,nsingle
cosph2(k)=0.0d0
sinph2(k)=0.0d0
phii1*rad2deg,ethetai
! lprn1=.false.
etheta=etheta+ethetai
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+ 'ebend',i,ethetai
if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
gloc(nphi+i-2,icg)=wang*dethetai
sumene1x,sumene2x,sumene3x,sumene4x,&
sumene1y,sumene2y,sumene3y,sumene4y,cossc,cossc1,&
cosfac2xx,sinfac2yy
-!el #define DEBUG
#ifdef DEBUG
real(kind=8) :: aincr,xxsave,sumenep,de_dxx_num,yysave,&
de_dyy_num,zzsave,de_dzz_num,costsave,sintsave,&
#ifdef DEBUG
write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i)
#endif
-!#undef DEBUG
!
!
cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
1 continue
enddo
-!el #undef DUBUG
return
end subroutine esc
!-----------------------------------------------------------------------------
etors_ii=0.0D0
if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 &
.or. itype(i).eq.ntyp1) cycle
- itori=itortyp(itype(i-2))
- itori1=itortyp(itype(i-1))
+ itori=itortyp(itype(i-2))
+ itori1=itortyp(itype(i-1))
phii=phi(i)
gloci=0.0D0
! Proline-Proline pair is a special case...
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
endif
- if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
'etor',i,etors_ii
if (lprn) &
write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') &
etors=0.0D0
do i=iphi_start,iphi_end
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
+ .or. itype(i-3).eq.ntyp1 &
.or. itype(i).eq.ntyp1) cycle
etors_ii=0.0D0
if (iabs(itype(i)).eq.20) then
! include 'COMMON.IOUNITS'
! include 'COMMON.FFIELD'
! include 'COMMON.TORCNSTR'
- real(kind=8) :: etors_d
+ real(kind=8) :: etors_d,etors_d_ii
logical :: lprn
!el local variables
integer :: i,j,k,l,itori,itori1,itori2,iblock
etors_d=0.0D0
! write(iout,*) "a tu??"
do i=iphid_start,iphid_end
+ etors_d_ii=0.0D0
if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 &
+ .or. itype(i-3).eq.ntyp1 &
.or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
itori=itortyp(itype(i-2))
itori1=itortyp(itype(i-1))
sinphi2=dsin(j*phii1)
etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+ &
v2cij*cosphi2+v2sij*sinphi2
+ if (energy_dec) etors_d_ii=etors_d_ii+ &
+ v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
enddo
sinphi1m2=dsin(l*phii-(k-l)*phii1)
etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+ &
v1sdij*sinphi1p2+v2sdij*sinphi1m2
+ if (energy_dec) etors_d_ii=etors_d_ii+ &
+ v1cdij*cosphi1p2+v2cdij*cosphi1m2+ &
+ v1sdij*sinphi1p2+v2sdij*sinphi1m2
gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2 &
-v1cdij*sinphi1p2-v2cdij*sinphi1m2)
gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2 &
-v1cdij*sinphi1p2+v2cdij*sinphi1m2)
enddo
enddo
+ if (energy_dec) write (iout,'(a6,i5,0pf7.3)') &
+ 'etor_d',i,etors_d_ii
gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
enddo
esccor_ii=0.0D0
isccori=isccortyp(itype(i-2))
isccori1=isccortyp(itype(i-1))
+
! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
phii=phi(i)
do intertyp=1,3 !intertyp
+ esccor_ii=0.0D0
!c Added 09 May 2012 (Adasko)
!c Intertyp means interaction type of backbone mainchain correlation:
! 1 = SC...Ca...Ca...Ca
v2ij=v2sccor(j,intertyp,isccori,isccori1)
cosphi=dcos(j*tauangle(intertyp,i))
sinphi=dsin(j*tauangle(intertyp,i))
+ if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+ if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)') &
+ 'esccor',i,intertyp,esccor_ii
! write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
if (lprn) &
#endif
#ifdef MPI
include 'mpif.h'
-!el#endif
+#endif
real(kind=8),dimension(3,nres) :: gradbufc,gradbufx,gradbufc_sum,&
gloc_scbuf !(3,maxres)
real(kind=8),dimension(4*nres) :: glocbuf !(4*maxres)
-#endif
+!#endif
!el local variables
integer :: i,j,k,ierror,ierr
real(kind=8) :: gvdwc_norm,gvdwc_scp_norm,gelc_norm,gvdwpp_norm,&
! " sigder",sigder
! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C print *,sss_ele_cut,'in sc_grad'
do k=1,3
dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
enddo
do k=1,3
- gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
+!C print *,'gg',k,gg(k)
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
gvdwx(k,i)=gvdwx(k,i)-gg(k) &
+(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
- +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv &
+ *sss_ele_cut
+
gvdwx(k,j)=gvdwx(k,j)+gg(k) &
+(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
- +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv &
+ *sss_ele_cut
+
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
! Check the gradient of the virtual-bond and SC vectors in the internal
! coordinates.
!
- aincr=1.0d-7
- aincr2=5.0d-8
+ aincr=1.0d-6
+ aincr2=5.0d-7
call cartder
write (iout,'(a)') '**************** dx/dalpha'
write (iout,'(a)')
nf=0
nfl=0
call zerograd
- aincr=1.0D-7
- print '(a)','CG processor',me,' calling CHECK_CART.'
+ aincr=1.0D-5
+ print '(a)','CG processor',me,' calling CHECK_CART.',aincr
nf=0
icall=0
call geom_to_var(nvar,x)
enddo
return
end subroutine check_ecart
+#ifdef CARGRAD
!-----------------------------------------------------------------------------
subroutine check_ecartint
! Check the gradient of the energy in Cartesian coordinates.
!el integer :: icall
!el common /srutu/ icall
real(kind=8),dimension(6) :: ggg,ggg1
- real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx,ddc1,ddcn
real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
- real(kind=8),dimension(3) :: dcnorm_safe,dxnorm_safe
+ real(kind=8),dimension(3) :: dcnorm_safe1,dcnorm_safe2,dxnorm_safe
real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
real(kind=8),dimension(0:n_ene) :: energia,energia1
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
+ write (iout,*) "Before geom_to_var"
call geom_to_var(nvar,x)
+ write (iout,*) "after geom_to_var"
+ write (iout,*) "split_ene ",split_ene
+ call flush(iout)
if (.not.split_ene) then
-write(iout,*) 'Calling CHECK_ECARTINT if'
+ write(iout,*) 'Calling CHECK_ECARTINT if'
call etotal(energia)
!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
etot=energia(0)
+ write (iout,*) "etot",etot
+ call flush(iout)
!el call enerprint(energia)
!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
call flush(iout)
enddo
endif
write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
+! do i=1,nres
+ do i=nnt,nct
+ do j=1,3
+ if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
+ if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
+ ddc(j)=c(j,i)
+ ddx(j)=c(j,i+nres)
+ dcnorm_safe1(j)=dc_norm(j,i-1)
+ dcnorm_safe2(j)=dc_norm(j,i)
+ dxnorm_safe(j)=dc_norm(j,i+nres)
+ enddo
+ do j=1,3
+ c(j,i)=ddc(j)+aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot1=energia1(0)
+ write (iout,*) "ij",i,j," etot1",etot1
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
+ c(j,i)=ddc(j)-aincr
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot2=energia1(0)
+ write (iout,*) "ij",i,j," etot2",etot2
+ ggg(j)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+! write (iout,*) "etot21",etot21," etot22",etot22
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i)=ddc(j)
+ if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
+ if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
+ if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i-1)=dcnorm_safe1(j)
+ dc_norm(j,i)=dcnorm_safe2(j)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ enddo
+ do j=1,3
+ c(j,i+nres)=ddx(j)+aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot1=energia1(0)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot11=energia1(0)
+ call etotal_short(energia1)
+ etot12=energia1(0)
+ endif
+!- end split gradient
+ c(j,i+nres)=ddx(j)-aincr
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ call int_from_cart1(.false.)
+ if (.not.split_ene) then
+ call etotal(energia1)
+ etot2=energia1(0)
+ ggg(j+3)=(etot1-etot2)/(2*aincr)
+ else
+!- split gradient
+ call etotal_long(energia1)
+ etot21=energia1(0)
+ ggg(j+3)=(etot11-etot21)/(2*aincr)
+ call etotal_short(energia1)
+ etot22=energia1(0)
+ ggg1(j+3)=(etot12-etot22)/(2*aincr)
+!- end split gradient
+ endif
+! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
+ c(j,i+nres)=ddx(j)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dxnorm_safe(j)
+ call int_from_cart1(.false.)
+ enddo
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
+ if (split_ene) then
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),&
+ k=1,6)
+ write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') &
+ i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),&
+ ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
+ endif
+ enddo
+ return
+ end subroutine check_ecartint
+#else
+!-----------------------------------------------------------------------------
+ subroutine check_ecartint
+! Check the gradient of the energy in Cartesian coordinates.
+ use io_base, only: intout
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.VAR'
+! include 'COMMON.CONTACTS'
+! include 'COMMON.MD'
+! include 'COMMON.LOCAL'
+! include 'COMMON.SPLITELE'
+ use comm_srutu
+!el integer :: icall
+!el common /srutu/ icall
+ real(kind=8),dimension(6) :: ggg,ggg1
+ real(kind=8),dimension(3) :: cc,xx,ddc,ddx
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(3) :: dcnorm_safe,dxnorm_safe
+ real(kind=8),dimension(6,0:nres) :: grad_s,grad_s1 !(6,0:maxres)
+ real(kind=8),dimension(nres) :: phi_temp,theta_temp,alph_temp,omeg_temp !(maxres)
+ real(kind=8),dimension(0:n_ene) :: energia,energia1
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+!EL external fdum
+ integer :: i,j,k,nf
+ real(kind=8) :: rlambd,aincr,etot,etot1,etot11,etot12,etot2,&
+ etot21,etot22
+ r_cut=2.0d0
+ rlambd=0.3d0
+ icg=1
+ nf=0
+ nfl=0
+ call intout
+! call intcartderiv
+! call checkintcartgrad
+ call zerograd
+ aincr=2.0D-5
+ write(iout,*) 'Calling CHECK_ECARTINT.',aincr
+ nf=0
+ icall=0
+ call geom_to_var(nvar,x)
+ if (.not.split_ene) then
+ call etotal(energia)
+ etot=energia(0)
+!el call enerprint(energia)
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ do i=1,nres
+ write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ else
+!- split gradient check
+ call zerograd
+ call etotal_long(energia)
+!el call enerprint(energia)
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ write (iout,*) "longrange grad"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s(j,i)=gcart(j,i)
+ grad_s(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ call zerograd
+ call etotal_short(energia)
+!el call enerprint(energia)
+ call flush(iout)
+ write (iout,*) "enter cartgrad"
+ call flush(iout)
+ call cartgrad
+ write (iout,*) "exit cartgrad"
+ call flush(iout)
+ icall =1
+ write (iout,*) "shortrange grad"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+ do j=1,3
+ grad_s1(j,0)=gcart(j,0)
+ enddo
+ do i=1,nres
+ do j=1,3
+ grad_s1(j,i)=gcart(j,i)
+ grad_s1(j+3,i)=gxcart(j,i)
+ enddo
+ enddo
+ endif
+ write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
do i=0,nres
do j=1,3
xx(j)=c(j,i+nres)
enddo
return
end subroutine check_ecartint
+#endif
!-----------------------------------------------------------------------------
subroutine check_eint
! Check the gradient of energy in internal coordinates.
call zerograd
aincr=1.0D-7
print '(a)','Calling CHECK_INT.'
-write(iout,*) 'Calling CHECK_INT.'
nf=0
nfl=0
icg=1
call geom_to_var(nvar,x)
call var_to_geom(nvar,x)
call chainbuild
-write(iout,*) 'Calling CHECK_INT.'
icall=1
print *,'ICG=',ICG
call etotal(energia)
nfl=3
!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
- write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
+!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
icall=1
do i=1,nvar
xi=x(i)
i,key,ii,gg(i),gana(i),&
100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
enddo
-write(iout,*) "jestesmy sobie w check eint!!"
return
end subroutine check_eint
!-----------------------------------------------------------------------------
endif
return
end function sscale
+ real(kind=8) function sscale_grad(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale_grad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscale_grad=0d0
+ endif
+ return
+ end function sscale_grad
+
+!!!!!!!!!! PBCSCALE
+ real(kind=8) function sscale_ele(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscale_ele=1.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_ele=0d0
+ endif
+ return
+ end function sscale_ele
+
+ real(kind=8) function sscagrad_ele(r)
+ real(kind=8) :: r,gamm
+! include "COMMON.SPLITELE"
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscagrad_ele=0.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+ else
+ sscagrad_ele=0.0d0
+ endif
+ return
+ end function sscagrad_ele
+!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------------
subroutine elj_long(evdw)
!
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
- real(kind=8) :: sss,e1,e2,evdw
+ real(kind=8) :: sss,e1,e2,evdw,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ 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
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+! Searching for nearest neighbour
+ 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
+
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
! Calculate angle-dependent terms of energy and contributions to their
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*(1.0d0-sss)
+ evdw=evdw+evdwij*(1.0d0-sss)*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_long"
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
- real(kind=8) :: sss,e1,e2,evdw,rij_shift
+ real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ 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
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+! xj=c(1,nres+j)-xi
+! yj=c(2,nres+j)-yi
+! zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+! Searching for nearest neighbour
+ 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
+
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ if (sss_ele_cut.le.0.0) cycle
if (sss.gt.0.0d0) then
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
+ evdw=evdw+evdwij*sss*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'evdw',i,j,evdwij
+! if (energy_dec) write (iout,*) &
+! 'evdw',i,j,evdwij,"egb_short"
! Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ /sigmaii(itypi,itypj))
+
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
! Contribution to the local-electrostatic energy coming from the i-j pair
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3) &
+a33*muij(4)
-!d write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
'eelloc',i,j,eel_loc_ij
-! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
+! write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) !d
eel_loc=eel_loc+eel_loc_ij
! Partial derivatives in virtual-bond dihedral angles gamma
dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
enddo
do k=1,3
- gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac&
+ *sss_ele_cut
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
gvdwx(k,i)=gvdwx(k,i)-gg(k) &
+(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
- +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac&
+ *sss_ele_cut
gvdwx(k,j)=gvdwx(k,j)+gg(k) &
+(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
- +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac&
+ *sss_ele_cut
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
! & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
!
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
- use MD_data, only: totT
+ use MD_data, only: totT,usampl,eq_time
#ifndef ISNAN
external proc_proc
#ifdef WINPGI
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
use energy_data
- use MD_data, only: totT
+ use MD_data, only: totT,usampl,eq_time
#ifdef MPI
include 'mpif.h'
#endif
call sum_gradient
#ifdef TIMING
#endif
+!el write (iout,*) "After sum_gradient"
#ifdef DEBUG
+!el write (iout,*) "After sum_gradient"
do i=1,nres-1
write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
enddo
endif
+!elwrite (iout,*) "After sum_gradient"
#ifdef TIMING
time01=MPI_Wtime()
#endif
call intcartderiv
+!elwrite (iout,*) "After sum_gradient"
#ifdef TIMING
time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
#endif
(gxcart(j,i),j=1,3)
enddo
#endif
+#ifdef CARGRAD
+#ifdef DEBUG
+ write (iout,*) "CARGRAD"
+#endif
+ do i=nres,1,-1
+ do j=1,3
+ gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+ enddo
+! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+ enddo
+! Correction: dummy residues
+ if (nnt.gt.1) then
+ do j=1,3
+! gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+ gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+ enddo
+ endif
+ if (nct.lt.nres) then
+ do j=1,3
+! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+ gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+ enddo
+ endif
+#endif
#ifdef TIMING
time_cartgrad=time_cartgrad+MPI_Wtime()-time00
#endif
! Obtaining the gamma derivatives from sine derivative
if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. &
phi(i).gt.pi34.and.phi(i).le.pi.or. &
- phi(i).gt.-pi.and.phi(i).le.-pi34) then
+ phi(i).ge.-pi.and.phi(i).le.-pi34) then
call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
cost=dcos(theta(i))
cost1=dcos(omicron(2,i-1))
cosg=dcos(tauangle(1,i))
+!elwrite(iout,*) " vecpr5",i,nres
do j=1,3
+!elwrite(iout,*) " vecpreee",i,nres,j,i-2+nres
+!elwrite(iout,*) " vecpr5",dc_norm2(1,1)
dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
enddo
integer :: kstart,kend,lstart,lend,idummy
real(kind=8) :: delta=1.0d-7
!el local variables
- integer :: i,ii,j
+ integer :: i,ii,j
! real(kind=8) ::
! For the backbone
do i=0,nres-1
real(kind=8) :: deps,ssx0,ljx0
!-------END TESTING CODE
+ eij=0.0d0
i=resi
j=resj
endif
if (havebond) then
-#ifndef CLUST
-#ifndef WHAM
+!#ifndef CLUST
+!#ifndef WHAM
! if (dyn_ssbond_ij(i,j).eq.1.0d300) then
! write(iout,'(a15,f12.2,f8.1,2i5)')
! & "SSBOND_E_FORM",totT,t_bath,i,j
! endif
-#endif
-#endif
+!#endif
+!#endif
dyn_ssbond_ij(i,j)=eij
else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
dyn_ssbond_ij(i,j)=1.0d300
-#ifndef CLUST
-#ifndef WHAM
+!#ifndef CLUST
+!#ifndef WHAM
! write(iout,'(a15,f12.2,f8.1,2i5)')
! & "SSBOND_E_BREAK",totT,t_bath,i,j
-#endif
-#endif
+!#endif
+!#endif
endif
!-------TESTING CODE
- if (checkstop) then
+!el if (checkstop) then
if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') &
"CHECKSTOP",rij,eij,ed
echeck(jcheck)=eij
- endif
+!el endif
enddo
if (checkstop) then
write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
! include 'COMMON.CHAIN'
! include 'COMMON.IOUNITS'
! include 'COMMON.SETUP'
-#ifndef CLUST
-#ifndef WHAM
! include 'COMMON.MD'
-#endif
-#endif
! Local variables
real(kind=8) :: emin
integer :: i,j,imin,ierr
!-----------------------------------------------------------------------------
#ifdef WHAM
subroutine read_ssHist
- implicit none
+! implicit none
! Includes
! include 'DIMENSIONS'
! include "DIMENSIONS.FREE"
!-----------------------------------------------------------------------------
subroutine alloc_ener_arrays
!EL Allocation of arrays used by module energy
-
+ use MD_data, only: mset
!el local variables
integer :: i,j
maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond
!----------------------
! arrays in subroutine init_int_table
+!el#ifdef MPI
+!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el#endif
allocate(nint_gr(nres))
allocate(nscp_gr(nres))
allocate(ielstart(nres))
- allocate(ielend(nres)) !(maxres)
+ allocate(ielend(nres))
+!(maxres)
allocate(istart(nres,maxint_gr))
- allocate(iend(nres,maxint_gr)) !(maxres,maxint_gr)
+ allocate(iend(nres,maxint_gr))
+!(maxres,maxint_gr)
allocate(iscpstart(nres,maxint_gr))
- allocate(iscpend(nres,maxint_gr)) !(maxres,maxint_gr)
+ allocate(iscpend(nres,maxint_gr))
+!(maxres,maxint_gr)
allocate(ielstart_vdw(nres))
- allocate(ielend_vdw(nres)) !(maxres)
+ allocate(ielend_vdw(nres))
+!(maxres)
- allocate(lentyp(0:nfgtasks-1)) !(0:maxprocs-1)
+ allocate(lentyp(0:nfgtasks-1))
+!(0:maxprocs-1)
!----------------------
! commom.contacts
! common /contacts/
if(.not.allocated(icont_ref)) allocate(icont_ref(2,maxcont))
- allocate(icont(2,maxcont)) !(2,maxcont)
+ allocate(icont(2,maxcont))
+!(2,maxcont)
! common /contacts1/
- allocate(num_cont(0:nres+4)) !(maxres)
- allocate(jcont(maxconts,nres)) !(maxconts,maxres)
- allocate(facont(maxconts,nres)) !(maxconts,maxres)
- allocate(gacont(3,maxconts,nres)) !(3,maxconts,maxres)
+ allocate(num_cont(0:nres+4))
+!(maxres)
+ allocate(jcont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(facont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(gacont(3,maxconts,nres))
+!(3,maxconts,maxres)
! common /contacts_hb/
allocate(gacontp_hb1(3,maxconts,nres))
allocate(gacontp_hb2(3,maxconts,nres))
allocate(gacontm_hb2(3,maxconts,nres))
allocate(gacontm_hb3(3,maxconts,nres))
allocate(gacont_hbr(3,maxconts,nres))
- allocate(grij_hb_cont(3,maxconts,nres)) !(3,maxconts,maxres)
+ allocate(grij_hb_cont(3,maxconts,nres))
+!(3,maxconts,maxres)
allocate(facont_hb(maxconts,nres))
allocate(ees0p(maxconts,nres))
allocate(ees0m(maxconts,nres))
- allocate(d_cont(maxconts,nres)) !(maxconts,maxres)
- allocate(num_cont_hb(nres)) !(maxres)
- allocate(jcont_hb(maxconts,nres)) !(maxconts,maxres)
+ allocate(d_cont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(num_cont_hb(nres))
+!(maxres)
+ allocate(jcont_hb(maxconts,nres))
+!(maxconts,maxres)
! common /rotat/
allocate(Ug(2,2,nres))
allocate(Ugder(2,2,nres))
allocate(Ug2(2,2,nres))
- allocate(Ug2der(2,2,nres)) !(2,2,maxres)
+ allocate(Ug2der(2,2,nres))
+!(2,2,maxres)
allocate(obrot(2,nres))
allocate(obrot2(2,nres))
allocate(obrot_der(2,nres))
- allocate(obrot2_der(2,nres)) !(2,maxres)
+ allocate(obrot2_der(2,nres))
+!(2,maxres)
! common /precomp1/
allocate(mu(2,nres))
allocate(muder(2,nres))
allocate(Ub2(2,nres))
+ Ub2(1,:)=0.0d0
+ Ub2(2,:)=0.0d0
allocate(Ub2der(2,nres))
allocate(Ctobr(2,nres))
allocate(Ctobrder(2,nres))
allocate(Dtobr2(2,nres))
- allocate(Dtobr2der(2,nres)) !(2,maxres)
+ allocate(Dtobr2der(2,nres))
+!(2,maxres)
allocate(EUg(2,2,nres))
allocate(EUgder(2,2,nres))
allocate(CUg(2,2,nres))
allocate(DUg(2,2,nres))
allocate(Dugder(2,2,nres))
allocate(DtUg2(2,2,nres))
- allocate(DtUg2der(2,2,nres)) !(2,2,maxres)
+ allocate(DtUg2der(2,2,nres))
+!(2,2,maxres)
! common /precomp2/
allocate(Ug2Db1t(2,nres))
allocate(Ug2Db1tder(2,nres))
allocate(CUgb2(2,nres))
- allocate(CUgb2der(2,nres)) !(2,maxres)
+ allocate(CUgb2der(2,nres))
+!(2,maxres)
allocate(EUgC(2,2,nres))
allocate(EUgCder(2,2,nres))
allocate(EUgD(2,2,nres))
allocate(EUgDder(2,2,nres))
allocate(DtUg2EUg(2,2,nres))
- allocate(Ug2DtEUg(2,2,nres)) !(2,2,maxres)
+ allocate(Ug2DtEUg(2,2,nres))
+!(2,2,maxres)
allocate(Ug2DtEUgder(2,2,2,nres))
- allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres)
+ allocate(DtUg2EUgder(2,2,2,nres))
+!(2,2,2,maxres)
! common /rotat_old/
allocate(costab(nres))
allocate(sintab(nres))
allocate(costab2(nres))
- allocate(sintab2(nres)) !(maxres)
+ allocate(sintab2(nres))
+!(maxres)
! common /dipmat/
allocate(a_chuj(2,2,maxconts,nres))
!(2,2,maxconts,maxres)(maxconts=maxres/4)
allocate(ncont_sent(nres))
allocate(ncont_recv(nres))
- allocate(iat_sent(nres)) !(maxres)
+ allocate(iat_sent(nres))
+!(maxres)
allocate(iint_sent(4,nres,nres))
- allocate(iint_sent_local(4,nres,nres)) !(4,maxres,maxres)
+ allocate(iint_sent_local(4,nres,nres))
+!(4,maxres,maxres)
allocate(iturn3_sent(4,0:nres+4))
allocate(iturn4_sent(4,0:nres+4))
allocate(iturn3_sent_local(4,nres))
- allocate(iturn4_sent_local(4,nres)) !(4,maxres)
+ allocate(iturn4_sent_local(4,nres))
+!(4,maxres)
allocate(itask_cont_from(0:nfgtasks-1))
- allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
+ allocate(itask_cont_to(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+
+
!----------------------
! commom.deriv;
! common /derivat/
allocate(dcdv(6,maxdim))
- allocate(dxdv(6,maxdim)) !(6,maxdim)
- allocate(dxds(6,nres)) !(6,maxres)
+ allocate(dxdv(6,maxdim))
+!(6,maxdim)
+ allocate(dxds(6,nres))
+!(6,maxres)
allocate(gradx(3,nres,0:2))
- allocate(gradc(3,nres,0:2)) !(3,maxres,2)
+ allocate(gradc(3,nres,0:2))
+!(3,maxres,2)
allocate(gvdwx(3,nres))
allocate(gvdwc(3,nres))
allocate(gelc(3,nres))
allocate(gcorr6_turn_long(3,nres))
allocate(gradxorr(3,nres))
allocate(gradcorr5(3,nres))
- allocate(gradcorr6(3,nres)) !(3,maxres)
+ allocate(gradcorr6(3,nres))
+!(3,maxres)
allocate(gloc(0:maxvar,0:2))
- allocate(gloc_x(0:maxvar,2)) !(maxvar,2)
+ allocate(gloc_x(0:maxvar,2))
+!(maxvar,2)
allocate(gel_loc(3,nres))
allocate(gel_loc_long(3,nres))
allocate(gcorr3_turn(3,nres))
allocate(gcorr4_turn(3,nres))
allocate(gcorr6_turn(3,nres))
allocate(gradb(3,nres))
- allocate(gradbx(3,nres)) !(3,maxres)
+ allocate(gradbx(3,nres))
+!(3,maxres)
allocate(gel_loc_loc(maxvar))
allocate(gel_loc_turn3(maxvar))
allocate(gel_loc_turn4(maxvar))
allocate(gel_loc_turn6(maxvar))
allocate(gcorr_loc(maxvar))
allocate(g_corr5_loc(maxvar))
- allocate(g_corr6_loc(maxvar)) !(maxvar)
+ allocate(g_corr6_loc(maxvar))
+!(maxvar)
allocate(gsccorc(3,nres))
- allocate(gsccorx(3,nres)) !(3,maxres)
- allocate(gsccor_loc(nres)) !(maxres)
- allocate(dtheta(3,2,nres)) !(3,2,maxres)
+ allocate(gsccorx(3,nres))
+!(3,maxres)
+ allocate(gsccor_loc(nres))
+!(maxres)
+ allocate(dtheta(3,2,nres))
+!(3,2,maxres)
allocate(gscloc(3,nres))
- allocate(gsclocx(3,nres)) !(3,maxres)
+ allocate(gsclocx(3,nres))
+!(3,maxres)
allocate(dphi(3,3,nres))
allocate(dalpha(3,3,nres))
- allocate(domega(3,3,nres)) !(3,3,maxres)
+ allocate(domega(3,3,nres))
+!(3,3,maxres)
! common /deriv_scloc/
allocate(dXX_C1tab(3,nres))
allocate(dYY_C1tab(3,nres))
allocate(dZZ_Ctab(3,nres))
allocate(dXX_XYZtab(3,nres))
allocate(dYY_XYZtab(3,nres))
- allocate(dZZ_XYZtab(3,nres)) !(3,maxres)
+ allocate(dZZ_XYZtab(3,nres))
+!(3,maxres)
! common /mpgrad/
allocate(jgrad_start(nres))
- allocate(jgrad_end(nres)) !(maxres)
+ allocate(jgrad_end(nres))
+!(maxres)
+!----------------------
! common /indices/
allocate(ibond_displ(0:nfgtasks-1))
allocate(iset_displ(0:nfgtasks-1))
allocate(iset_count(0:nfgtasks-1))
allocate(iint_count(0:nfgtasks-1))
- allocate(iint_displ(0:nfgtasks-1)) !(0:max_fg_procs-1)
+ allocate(iint_displ(0:nfgtasks-1))
+!(0:max_fg_procs-1)
!----------------------
! common.MD
! common /mdgrad/
allocate(gcart(3,0:nres))
- allocate(gxcart(3,0:nres)) !(3,0:MAXRES)
+ allocate(gxcart(3,0:nres))
+!(3,0:MAXRES)
allocate(gradcag(3,nres))
- allocate(gradxag(3,nres)) !(3,MAXRES)
+ allocate(gradxag(3,nres))
+!(3,MAXRES)
! common /back_constr/
!el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back)
allocate(dutheta(nres))
- allocate(dugamma(nres)) !(maxres)
+ allocate(dugamma(nres))
+!(maxres)
allocate(duscdiff(3,nres))
- allocate(duscdiffx(3,nres)) !(3,maxres)
+ allocate(duscdiffx(3,nres))
+!(3,maxres)
!el i io:read_fragments
! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)
! allocate((:,:,:),allocatable :: ifrag_back !(3,maxfrag_back,maxprocs/20)
! allocate(qinfrag(50,nprocs/20),wfrag(50,nprocs/20)) !(50,maxprocs/20)
! allocate(qinpair(100,nprocs/20),wpair(100,nprocs/20)) !(100,maxprocs/20)
allocate(mset(0:nprocs)) !(maxprocs/20)
- do i=0,nprocs
- mset(i)=0
- enddo
+ mset(:)=0
! allocate(ifrag(2,50,nprocs/20)) !(2,50,maxprocs/20)
! allocate(ipair(2,100,nprocs/20)) !(2,100,maxprocs/20)
allocate(dUdconst(3,0:nres))
allocate(dUdxconst(3,0:nres))
allocate(dqwol(3,0:nres))
- allocate(dxqwol(3,0:nres)) !(3,0:MAXRES)
+ allocate(dxqwol(3,0:nres))
+!(3,0:MAXRES)
!----------------------
! common.sbridge
! common /sbridge/ in io_common: read_bridge
!el integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
! common /dyn_ssbond/
! and side-chain vectors in theta or phi.
- allocate(dyn_ssbond_ij(0:nres+4,0:nres+4)) !(maxres,maxres)
- do i=1,nres
- do j=i+1,nres
- dyn_ssbond_ij(i,j)=1.0d300
- enddo
- enddo
+ allocate(dyn_ssbond_ij(0:nres+4,0:nres+4))
+!(maxres,maxres)
+! do i=1,nres
+! do j=i+1,nres
+ dyn_ssbond_ij(:,:)=1.0d300
+! enddo
+! enddo
if (nss.gt.0) then
- allocate(idssb(nss),jdssb(nss)) !(maxdim)
+ allocate(idssb(nss),jdssb(nss))
+!(maxdim)
endif
- allocate(dyn_ss_mask(nres)) !(maxres)
- do i=1,nres
- dyn_ss_mask(i)=.false.
- enddo
+ allocate(dyn_ss_mask(nres))
+!(maxres)
+ dyn_ss_mask(:)=.false.
!----------------------
! common.sccor
! Parameters of the SCCOR term
! allocate(vlor2sccor(maxterm_sccor,20,20))
! allocate(vlor3sccor(maxterm_sccor,20,20)) !(maxterm_sccor,20,20)
!----------------
- allocate(gloc_sc(3,0:2*nres,0:10)) !(3,0:maxres2,10)maxres2=2*maxres
+ allocate(gloc_sc(3,0:2*nres,0:10))
+!(3,0:maxres2,10)maxres2=2*maxres
allocate(dcostau(3,3,3,2*nres))
allocate(dsintau(3,3,3,2*nres))
allocate(dtauangle(3,3,3,2*nres))
allocate(dcosomicron(3,3,3,2*nres))
- allocate(domicron(3,3,3,2*nres)) !(3,3,3,maxres2)maxres2=2*maxres
-!----------------------
-! common.scrot
-! Parameters of the SC rotamers (local) term
-! common/scrot/ in io_conf: parmread
-! allocate((:,:),allocatable :: sc_parmin !(maxsccoef,ntyp)
-!----------------------
-! common.torcnstr
-! common /torcnstr/
-!el in io_conf:molread
-! allocate((:),allocatable :: idih_constr,idih_nconstr !(maxdih_constr)
-! allocate((:),allocatable :: phi0,drange !(maxdih_constr)
-!----------------------
-! common.torsion
-! common/torsion/ in io_conf: parmread
-! allocate((:,:,:),allocatable :: v0 !(-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:),allocatable :: v1,v2 !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:),allocatable :: vlor1 !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
-! allocate((:,:,:),allocatable :: vlor2,vlor3 !(maxlor,maxtor,maxtor)
-! allocate((:),allocatable :: itortyp !(-ntyp1:ntyp1)
-! allocate((:,:,:),allocatable :: nterm,nlor !(-maxtor:maxtor,-maxtor:maxtor,2)
-!
-! common /torsiond/ in io_conf: parmread
-! allocate((:,:,:,:,:,:),allocatable :: v1c,v1s
- !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:,:,:),allocatable :: v2c,v2s
- !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! allocate((:,:,:,:),allocatable :: ntermd_1,ntermd_2
- !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
-! common/fourier/ in io_conf: parmread
-! allocate((:,:),allocatable :: b1,b2,&
-! b1tilde !(2,-maxtor:maxtor)
-! allocate((:,:,:),allocatable :: cc,dd,ee,&
-! ctilde,dtilde !(2,2,-maxtor:maxtor)
+ allocate(domicron(3,3,3,2*nres))
+!(3,3,3,maxres2)maxres2=2*maxres
!----------------------
! common.var
! common /restr/
- allocate(varall(maxvar)) !(maxvar)(maxvar=6*maxres)
+ allocate(varall(maxvar))
+!(maxvar)(maxvar=6*maxres)
allocate(mask_theta(nres))
allocate(mask_phi(nres))
- allocate(mask_side(nres)) !(maxres)
+ allocate(mask_side(nres))
+!(maxres)
!----------------------
! common.vectors
! common /vectors/
allocate(uy(3,nres))
- allocate(uz(3,nres)) !(3,maxres)
+ allocate(uz(3,nres))
+!(3,maxres)
allocate(uygrad(3,3,2,nres))
- allocate(uzgrad(3,3,2,nres)) !(3,3,2,maxres)
+ allocate(uzgrad(3,3,2,nres))
+!(3,3,2,maxres)
return
end subroutine alloc_ener_arrays
projects / unres4.git / blobdiff