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.ntyp1) 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)) if (itypi.gt.ntyp .or. itypi.le.0) cycle itypi1=iabs(itype(i+1)) 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 if ((zi.gt.bordlipbot) &.and.(zi.lt.bordliptop)) then C the energy transfer exist if (zi.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipi=sscalelip(fracinbuf) ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick elseif (zi.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) sslipi=sscalelip(fracinbuf) ssgradlipi=sscagradlip(fracinbuf)/lipbufthick else sslipi=1.0d0 ssgradlipi=0.0 endif else sslipi=0.0d0 ssgradlipi=0.0 endif 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)) if (itypj.gt.ntyp .or. itypj.le.0) cycle 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 C xj=c(1,nres+j)-xi C yj=c(2,nres+j)-yi C zj=c(3,nres+j)-zi xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+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 if ((zj.gt.bordlipbot) &.and.(zj.lt.bordliptop)) then C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- & ((positi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick elseif (zi.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else sslipj=1.0d0 ssgradlipj=0.0 endif else sslipj=0.0d0 ssgradlipj=0.0 endif aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 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)/sigma(itypi,itypj)) sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj)) 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 e2=fac*bb 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/bb)**(1.0D0/6.0D0) epsi=bb**2/aa 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 fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij C Calculate the radial part of the gradient gg(1)=xj*fac gg(2)=yj*fac gg(3)=zj*fac gg_lipi(3)=ssgradlipi*evdwij gg_lipj(3)=ssgradlipj*evdwij C Calculate angular part of the gradient. call sc_grad ENDIF enddo ! j enddo ! iint enddo ! i end C-----------------------------------------------------------------------------