+++ /dev/null
- 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-----------------------------------------------------------------------------
projects / unres.git / blobdiff