X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Fsc_move.F;fp=source%2Funres%2Fsrc_CSA_DiL%2Fsc_move.F;h=0000000000000000000000000000000000000000;hp=74e9bf2d00a941d45451cd8ba21f24b9f4dbb78b;hb=de4bc5453ea46e111d936cb85e1758ed21c08fcd;hpb=b75425747e3e2b448ca5e0ef8367712e1f339124 diff --git a/source/unres/src_CSA_DiL/sc_move.F b/source/unres/src_CSA_DiL/sc_move.F deleted file mode 100644 index 74e9bf2..0000000 --- a/source/unres/src_CSA_DiL/sc_move.F +++ /dev/null @@ -1,823 +0,0 @@ - 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(iabs(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=iabs(itype(i)) - itypi1=iabs(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=iabs(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-----------------------------------------------------------------------------