--- /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'
+ include 'mpif.h'
+ 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))
+ 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