X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Fsc_move.F;fp=source%2Funres%2Fsrc_MD-NEWSC%2Fsc_move.F;h=b6837fde271a4f35fe8a920d4432dcfdb3d3b80b;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-NEWSC/sc_move.F b/source/unres/src_MD-NEWSC/sc_move.F new file mode 100644 index 0000000..b6837fd --- /dev/null +++ b/source/unres/src_MD-NEWSC/sc_move.F @@ -0,0 +1,823 @@ + subroutine sc_move(n_start,n_end,n_maxtry,e_drop, + + n_fun,etot) +c Perform a quick search over side-chain arrangments (over +c residues n_start to n_end) for a given (frozen) CA trace +c Only side-chains are minimized (at most n_maxtry times each), +c not CA positions +c Stops if energy drops by e_drop, otherwise tries all residues +c in the given range +c If there is an energy drop, full minimization may be useful +c n_start, n_end CAN be modified by this routine, but only if +c out of bounds (n_start <= 1, n_end >= nres, n_start < n_end) +c NOTE: this move should never increase the energy +crc implicit none + +c Includes + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.HEADER' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.FFIELD' + +c External functions + integer iran_num + external iran_num + +c Input arguments + integer n_start,n_end,n_maxtry + double precision e_drop + +c Output arguments + integer n_fun + double precision etot + +c Local variables + double precision energy(0:n_ene) + double precision cur_alph(2:nres-1),cur_omeg(2:nres-1) + double precision orig_e,cur_e + integer n,n_steps,n_first,n_cur,n_tot,i + double precision orig_w(n_ene) + double precision wtime + + +c Set non side-chain weights to zero (minimization is faster) +c NOTE: e(2) does not actually depend on the side-chain, only CA + orig_w(2)=wscp + orig_w(3)=welec + orig_w(4)=wcorr + orig_w(5)=wcorr5 + orig_w(6)=wcorr6 + orig_w(7)=wel_loc + orig_w(8)=wturn3 + orig_w(9)=wturn4 + orig_w(10)=wturn6 + orig_w(11)=wang + orig_w(13)=wtor + orig_w(14)=wtor_d + orig_w(15)=wvdwpp + + wscp=0.D0 + welec=0.D0 + wcorr=0.D0 + wcorr5=0.D0 + wcorr6=0.D0 + wel_loc=0.D0 + wturn3=0.D0 + wturn4=0.D0 + wturn6=0.D0 + wang=0.D0 + wtor=0.D0 + wtor_d=0.D0 + wvdwpp=0.D0 + +c Make sure n_start, n_end are within proper range + if (n_start.lt.2) n_start=2 + if (n_end.gt.nres-1) n_end=nres-1 +crc if (n_start.lt.n_end) then + if (n_start.gt.n_end) then + n_start=2 + n_end=nres-1 + endif + +c Save the initial values of energy and coordinates +cd call chainbuild +cd call etotal(energy) +cd write (iout,*) 'start sc ene',energy(0) +cd call enerprint(energy(0)) +crc etot=energy(0) + n_fun=0 +crc orig_e=etot +crc cur_e=orig_e +crc do i=2,nres-1 +crc cur_alph(i)=alph(i) +crc cur_omeg(i)=omeg(i) +crc enddo + +ct wtime=MPI_WTIME() +c Try (one by one) all specified residues, starting from a +c random position in sequence +c Stop early if the energy has decreased by at least e_drop + n_tot=n_end-n_start+1 + n_first=iran_num(0,n_tot-1) + n_steps=0 + n=0 +crc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop) + do while (n.lt.n_tot) + n_cur=n_start+mod(n_first+n,n_tot) + call single_sc_move(n_cur,n_maxtry,e_drop, + + n_steps,n_fun,etot) +c If a lower energy was found, update the current structure... +crc if (etot.lt.cur_e) then +crc cur_e=etot +crc do i=2,nres-1 +crc cur_alph(i)=alph(i) +crc cur_omeg(i)=omeg(i) +crc enddo +crc else +c ...else revert to the previous one +crc etot=cur_e +crc do i=2,nres-1 +crc alph(i)=cur_alph(i) +crc omeg(i)=cur_omeg(i) +crc enddo +crc endif + n=n+1 +cd +cd call chainbuild +cd call etotal(energy) +cd print *,'running',n,energy(0) + enddo + +cd call chainbuild +cd call etotal(energy) +cd write (iout,*) 'end sc ene',energy(0) + +c Put the original weights back to calculate the full energy + wscp=orig_w(2) + welec=orig_w(3) + wcorr=orig_w(4) + wcorr5=orig_w(5) + wcorr6=orig_w(6) + wel_loc=orig_w(7) + wturn3=orig_w(8) + wturn4=orig_w(9) + wturn6=orig_w(10) + wang=orig_w(11) + wtor=orig_w(13) + wtor_d=orig_w(14) + wvdwpp=orig_w(15) + +crc n_fun=n_fun+1 +ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime + return + end + +c------------------------------------------------------------- + + subroutine single_sc_move(res_pick,n_maxtry,e_drop, + + n_steps,n_fun,e_sc) +c Perturb one side-chain (res_pick) and minimize the +c neighbouring region, keeping all CA's and non-neighbouring +c side-chains fixed +c Try until e_drop energy improvement is achieved, or n_maxtry +c attempts have been made +c At the start, e_sc should contain the side-chain-only energy(0) +c nsteps and nfun for this move are ADDED to n_steps and n_fun +crc implicit none + +c Includes + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.CHAIN' + include 'COMMON.MINIM' + include 'COMMON.FFIELD' + include 'COMMON.IOUNITS' + +c External functions + double precision dist + external dist + +c Input arguments + integer res_pick,n_maxtry + double precision e_drop + +c Input/Output arguments + integer n_steps,n_fun + double precision e_sc + +c Local variables + logical fail + integer i,j + integer nres_moved + integer iretcode,loc_nfun,orig_maxfun,n_try + double precision sc_dist,sc_dist_cutoff + double precision energy(0:n_ene),orig_e,cur_e + double precision evdw,escloc + double precision cur_alph(2:nres-1),cur_omeg(2:nres-1) + double precision var(maxvar) + + double precision orig_theta(1:nres),orig_phi(1:nres), + + orig_alph(1:nres),orig_omeg(1:nres) + + +c Define what is meant by "neighbouring side-chain" + sc_dist_cutoff=5.0D0 + +c Don't do glycine or ends + i=itype(res_pick) + if (i.eq.10 .or. i.eq.21) return + +c Freeze everything (later will relax only selected side-chains) + mask_r=.true. + do i=1,nres + mask_phi(i)=0 + mask_theta(i)=0 + mask_side(i)=0 + enddo + +c Find the neighbours of the side-chain to move +c and save initial variables +crc orig_e=e_sc +crc cur_e=orig_e + nres_moved=0 + do i=2,nres-1 +c Don't do glycine (itype(j)==10) + if (itype(i).ne.10) then + sc_dist=dist(nres+i,nres+res_pick) + else + sc_dist=sc_dist_cutoff + endif + if (sc_dist.lt.sc_dist_cutoff) then + nres_moved=nres_moved+1 + mask_side(i)=1 + cur_alph(i)=alph(i) + cur_omeg(i)=omeg(i) + endif + enddo + + call chainbuild + call egb1(evdw) + call esc(escloc) + e_sc=wsc*evdw+wscloc*escloc +cd call etotal(energy) +cd print *,'new ',(energy(k),k=0,n_ene) + orig_e=e_sc + cur_e=orig_e + + n_try=0 + do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop) +c Move the selected residue (don't worry if it fails) + call gen_side(itype(res_pick),theta(res_pick+1), + + alph(res_pick),omeg(res_pick),fail) + +c Minimize the side-chains starting from the new arrangement + call geom_to_var(nvar,var) + orig_maxfun=maxfun + maxfun=7 + +crc do i=1,nres +crc orig_theta(i)=theta(i) +crc orig_phi(i)=phi(i) +crc orig_alph(i)=alph(i) +crc orig_omeg(i)=omeg(i) +crc enddo + + call minimize_sc1(e_sc,var,iretcode,loc_nfun) + +cv write(*,'(2i3,2f12.5,2i3)') +cv & res_pick,nres_moved,orig_e,e_sc-cur_e, +cv & iretcode,loc_nfun + +c$$$ if (iretcode.eq.8) then +c$$$ write(iout,*)'Coordinates just after code 8' +c$$$ call chainbuild +c$$$ call all_varout +c$$$ call flush(iout) +c$$$ do i=1,nres +c$$$ theta(i)=orig_theta(i) +c$$$ phi(i)=orig_phi(i) +c$$$ alph(i)=orig_alph(i) +c$$$ omeg(i)=orig_omeg(i) +c$$$ enddo +c$$$ write(iout,*)'Coordinates just before code 8' +c$$$ call chainbuild +c$$$ call all_varout +c$$$ call flush(iout) +c$$$ endif + + n_fun=n_fun+loc_nfun + maxfun=orig_maxfun + call var_to_geom(nvar,var) + +c If a lower energy was found, update the current structure... + if (e_sc.lt.cur_e) then +cv call chainbuild +cv call etotal(energy) +cd call egb1(evdw) +cd call esc(escloc) +cd e_sc1=wsc*evdw+wscloc*escloc +cd print *,' new',e_sc1,energy(0) +cv print *,'new ',energy(0) +cd call enerprint(energy(0)) + cur_e=e_sc + do i=2,nres-1 + if (mask_side(i).eq.1) then + cur_alph(i)=alph(i) + cur_omeg(i)=omeg(i) + endif + enddo + else +c ...else revert to the previous one + e_sc=cur_e + do i=2,nres-1 + if (mask_side(i).eq.1) then + alph(i)=cur_alph(i) + omeg(i)=cur_omeg(i) + endif + enddo + endif + n_try=n_try+1 + + enddo + n_steps=n_steps+n_try + +c Reset the minimization mask_r to false + mask_r=.false. + + return + end + +c------------------------------------------------------------- + + subroutine sc_minimize(etot,iretcode,nfun) +c Minimizes side-chains only, leaving backbone frozen +crc implicit none + +c Includes + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.FFIELD' + +c Output arguments + double precision etot + integer iretcode,nfun + +c Local variables + integer i + double precision orig_w(n_ene),energy(0:n_ene) + double precision var(maxvar) + + +c Set non side-chain weights to zero (minimization is faster) +c NOTE: e(2) does not actually depend on the side-chain, only CA + orig_w(2)=wscp + orig_w(3)=welec + orig_w(4)=wcorr + orig_w(5)=wcorr5 + orig_w(6)=wcorr6 + orig_w(7)=wel_loc + orig_w(8)=wturn3 + orig_w(9)=wturn4 + orig_w(10)=wturn6 + orig_w(11)=wang + orig_w(13)=wtor + orig_w(14)=wtor_d + + wscp=0.D0 + welec=0.D0 + wcorr=0.D0 + wcorr5=0.D0 + wcorr6=0.D0 + wel_loc=0.D0 + wturn3=0.D0 + wturn4=0.D0 + wturn6=0.D0 + wang=0.D0 + wtor=0.D0 + wtor_d=0.D0 + +c Prepare to freeze backbone + do i=1,nres + mask_phi(i)=0 + mask_theta(i)=0 + mask_side(i)=1 + enddo + +c Minimize the side-chains + mask_r=.true. + call geom_to_var(nvar,var) + call minimize(etot,var,iretcode,nfun) + call var_to_geom(nvar,var) + mask_r=.false. + +c Put the original weights back and calculate the full energy + wscp=orig_w(2) + welec=orig_w(3) + wcorr=orig_w(4) + wcorr5=orig_w(5) + wcorr6=orig_w(6) + wel_loc=orig_w(7) + wturn3=orig_w(8) + wturn4=orig_w(9) + wturn6=orig_w(10) + wang=orig_w(11) + wtor=orig_w(13) + wtor_d=orig_w(14) + + call chainbuild + call etotal(energy) + etot=energy(0) + + return + end + +c------------------------------------------------------------- + subroutine minimize_sc1(etot,x,iretcode,nfun) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) + include 'COMMON.IOUNITS' + include 'COMMON.VAR' + include 'COMMON.GEO' + include 'COMMON.MINIM' + common /srutu/ icall + dimension iv(liv) + double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar) + double precision energia(0:n_ene) + external func,gradient,fdum + external func_restr1,grad_restr1 + logical not_done,change,reduce + common /przechowalnia/ v + + call deflt(2,iv,liv,lv,v) +* 12 means fresh start, dont call deflt + iv(1)=12 +* max num of fun calls + if (maxfun.eq.0) maxfun=500 + iv(17)=maxfun +* max num of iterations + if (maxmin.eq.0) maxmin=1000 + iv(18)=maxmin +* controls output + iv(19)=2 +* selects output unit +c iv(21)=iout + iv(21)=0 +* 1 means to print out result + iv(22)=0 +* 1 means to print out summary stats + iv(23)=0 +* 1 means to print initial x and d + iv(24)=0 +* min val for v(radfac) default is 0.1 + v(24)=0.1D0 +* max val for v(radfac) default is 4.0 + v(25)=2.0D0 +c v(25)=4.0D0 +* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease) +* the sumsl default is 0.1 + v(26)=0.1D0 +* false conv if (act fnctn decrease) .lt. v(34) +* the sumsl default is 100*machep + v(34)=v(34)/100.0D0 +* absolute convergence + if (tolf.eq.0.0D0) tolf=1.0D-4 + v(31)=tolf +* relative convergence + if (rtolf.eq.0.0D0) rtolf=1.0D-4 + v(32)=rtolf +* controls initial step size + v(35)=1.0D-1 +* large vals of d correspond to small components of step + do i=1,nphi + d(i)=1.0D-1 + enddo + do i=nphi+1,nvar + d(i)=1.0D-1 + enddo + IF (mask_r) THEN + call x2xx(x,xx,nvar_restr) + call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1, + & iv,liv,lv,v,idum,rdum,fdum) + call xx2x(x,xx) + ELSE + call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum) + ENDIF + etot=v(10) + iretcode=iv(1) + nfun=iv(6) + + return + end +************************************************************************ + subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.DERIV' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.FFIELD' + include 'COMMON.INTERACT' + include 'COMMON.TIME1' + common /chuju/ jjj + double precision energia(0:n_ene),evdw,escloc + integer jjj + double precision ufparm,e1,e2 + external ufparm + integer uiparm(1) + real*8 urparm(1) + dimension x(maxvar) + nfl=nf + icg=mod(nf,2)+1 + +#ifdef OSF +c Intercept NaNs in the coordinates, before calling etotal + x_sum=0.D0 + do i=1,n + x_sum=x_sum+x(i) + enddo + FOUND_NAN=.false. + if (x_sum.ne.x_sum) then + write(iout,*)" *** func_restr1 : Found NaN in coordinates" + f=1.0D+73 + FOUND_NAN=.true. + return + endif +#endif + + call var_to_geom_restr(n,x) + call zerograd + call chainbuild +cd write (iout,*) 'ETOTAL called from FUNC' + call egb1(evdw) + call esc(escloc) + f=wsc*evdw+wscloc*escloc +cd call etotal(energia(0)) +cd f=wsc*energia(1)+wscloc*energia(12) +cd print *,f,evdw,escloc,energia(0) +C +C Sum up the components of the Cartesian gradient. +C + do i=1,nct + do j=1,3 + gradx(j,i,icg)=wsc*gvdwx(j,i) + enddo + enddo + + return + end +c------------------------------------------------------- + subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.VAR' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + include 'COMMON.IOUNITS' + external ufparm + integer uiparm(1) + double precision urparm(1) + dimension x(maxvar),g(maxvar) + + icg=mod(nf,2)+1 + if (nf-nfl+1) 20,30,40 + 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm) +c write (iout,*) 'grad 20' + if (nf.eq.0) return + goto 40 + 30 call var_to_geom_restr(n,x) + call chainbuild +C +C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. +C + 40 call cartder +C +C Convert the Cartesian gradient into internal-coordinate gradient. +C + + ig=0 + ind=nres-2 + do i=2,nres-2 + IF (mask_phi(i+2).eq.1) THEN + gphii=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) + gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) + enddo + enddo + ig=ig+1 + g(ig)=gphii + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + + ind=0 + do i=1,nres-2 + IF (mask_theta(i+2).eq.1) THEN + ig=ig+1 + gthetai=0.0D0 + do j=i+1,nres-1 + ind=ind+1 + do k=1,3 + gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) + gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) + enddo + enddo + g(ig)=gthetai + ELSE + ind=ind+nres-1-i + ENDIF + enddo + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + galphai=0.0D0 + do k=1,3 + galphai=galphai+dxds(k,i)*gradx(k,i,icg) + enddo + g(ig)=galphai + ENDIF + endif + enddo + + + do i=2,nres-1 + if (itype(i).ne.10) then + IF (mask_side(i).eq.1) THEN + ig=ig+1 + gomegai=0.0D0 + do k=1,3 + gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) + enddo + g(ig)=gomegai + ENDIF + endif + enddo + +C +C Add the components corresponding to local energy terms. +C + + ig=0 + igall=0 + do i=4,nres + igall=igall+1 + if (mask_phi(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do i=3,nres + igall=igall+1 + if (mask_theta(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + enddo + + do ij=1,2 + do i=2,nres-1 + if (itype(i).ne.10) then + igall=igall+1 + if (mask_side(i).eq.1) then + ig=ig+1 + g(ig)=g(ig)+gloc(igall,icg) + endif + endif + enddo + enddo + +cd do i=1,ig +cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) +cd enddo + return + end +C----------------------------------------------------------------------------- + subroutine egb1(evdw) +C +C This subroutine calculates the interaction energy of nonbonded side chains +C assuming the Gay-Berne potential of interaction. +C + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + logical lprn + evdw=0.0D0 +c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon + evdw=0.0D0 + lprn=.false. +c if (icall.eq.0) lprn=.true. + ind=0 + do i=iatsc_s,iatsc_e + + + itypi=itype(i) + itypi1=itype(i+1) + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + dxi=dc_norm(1,nres+i) + dyi=dc_norm(2,nres+i) + dzi=dc_norm(3,nres+i) + 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) + IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN + ind=ind+1 + itypj=itype(j) + dscj_inv=dsc_inv(itypj) + sig0ij=sigma(itypi,itypj) + chi1=chi(itypi,itypj) + chi2=chi(itypj,itypi) + chi12=chi1*chi2 + chip1=chip(itypi) + chip2=chip(itypj) + chip12=chip1*chip2 + alf1=alp(itypi) + alf2=alp(itypj) + alf12=0.5D0*(alf1+alf2) +C For diagnostics only!!! +c chi1=0.0D0 +c chi2=0.0D0 +c chi12=0.0D0 +c chip1=0.0D0 +c chip2=0.0D0 +c chip12=0.0D0 +c alf1=0.0D0 +c alf2=0.0D0 +c alf12=0.0D0 + xj=c(1,nres+j)-xi + yj=c(2,nres+j)-yi + zj=c(3,nres+j)-zi + 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) +C Calculate angle-dependent terms of energy and contributions to their +C derivatives. + call sc_angular + sigsq=1.0D0/sigsq + sig=sig0ij*dsqrt(sigsq) + rij_shift=1.0D0/rij-sig+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 +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) + return + endif + sigder=-sig*sigsq +c--------------------------------------------------------------- + rij_shift=1.0D0/rij_shift + fac=rij_shift**expon + e1=fac*fac*aa(itypi,itypj) + e2=fac*bb(itypi,itypj) + evdwij=eps1*eps2rt*eps3rt*(e1+e2) + eps2der=evdwij*eps3rt + eps3der=evdwij*eps2rt + evdwij=evdwij*eps2rt*eps3rt + evdw=evdw+evdwij + if (lprn) then + sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) + epsi=bb(itypi,itypj)**2/aa(itypi,itypj) +cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') +cd & restyp(itypi),i,restyp(itypj),j, +cd & epsi,sigm,chi1,chi2,chip1,chip2, +cd & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, +cd & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, +cd & evdwij + endif + + if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') + & 'evdw',i,j,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 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. + call sc_grad + ENDIF + enddo ! j + enddo ! iint + enddo ! i + end +C-----------------------------------------------------------------------------