1 subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
3 c Perform a quick search over side-chain arrangments (over
4 c residues n_start to n_end) for a given (frozen) CA trace
5 c Only side-chains are minimized (at most n_maxtry times each),
7 c Stops if energy drops by e_drop, otherwise tries all residues
9 c If there is an energy drop, full minimization may be useful
10 c n_start, n_end CAN be modified by this routine, but only if
11 c out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
12 c NOTE: this move should never increase the energy
16 implicit real*8 (a-h,o-z)
23 include 'COMMON.HEADER'
24 include 'COMMON.IOUNITS'
25 include 'COMMON.CHAIN'
26 include 'COMMON.FFIELD'
33 integer n_start,n_end,n_maxtry
34 double precision e_drop
41 double precision energy(0:n_ene)
42 double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
43 double precision orig_e,cur_e
44 integer n,n_steps,n_first,n_cur,n_tot,i
45 double precision orig_w(n_ene)
46 double precision wtime
49 c Set non side-chain weights to zero (minimization is faster)
50 c NOTE: e(2) does not actually depend on the side-chain, only CA
79 c Make sure n_start, n_end are within proper range
80 if (n_start.lt.2) n_start=2
81 if (n_end.gt.nres-1) n_end=nres-1
82 crc if (n_start.lt.n_end) then
83 if (n_start.gt.n_end) then
88 c Save the initial values of energy and coordinates
90 cd call etotal(energy)
91 cd write (iout,*) 'start sc ene',energy(0)
92 cd call enerprint(energy(0))
98 crc cur_alph(i)=alph(i)
99 crc cur_omeg(i)=omeg(i)
103 c Try (one by one) all specified residues, starting from a
104 c random position in sequence
105 c Stop early if the energy has decreased by at least e_drop
106 n_tot=n_end-n_start+1
107 n_first=iran_num(0,n_tot-1)
110 crc do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
111 do while (n.lt.n_tot)
112 n_cur=n_start+mod(n_first+n,n_tot)
113 call single_sc_move(n_cur,n_maxtry,e_drop,
114 + n_steps,n_fun,etot)
115 c If a lower energy was found, update the current structure...
116 crc if (etot.lt.cur_e) then
119 crc cur_alph(i)=alph(i)
120 crc cur_omeg(i)=omeg(i)
123 c ...else revert to the previous one
126 crc alph(i)=cur_alph(i)
127 crc omeg(i)=cur_omeg(i)
133 cd call etotal(energy)
134 cd print *,'running',n,energy(0)
138 cd call etotal(energy)
139 cd write (iout,*) 'end sc ene',energy(0)
141 c Put the original weights back to calculate the full energy
157 ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
161 c-------------------------------------------------------------
163 subroutine single_sc_move(res_pick,n_maxtry,e_drop,
164 + n_steps,n_fun,e_sc)
165 c Perturb one side-chain (res_pick) and minimize the
166 c neighbouring region, keeping all CA's and non-neighbouring
168 c Try until e_drop energy improvement is achieved, or n_maxtry
169 c attempts have been made
170 c At the start, e_sc should contain the side-chain-only energy(0)
171 c nsteps and nfun for this move are ADDED to n_steps and n_fun
175 implicit real*8 (a-h,o-z)
178 include 'COMMON.INTERACT'
179 include 'COMMON.CHAIN'
180 include 'COMMON.MINIM'
181 include 'COMMON.FFIELD'
182 include 'COMMON.IOUNITS'
185 double precision dist
189 integer res_pick,n_maxtry
190 double precision e_drop
192 c Input/Output arguments
193 integer n_steps,n_fun
194 double precision e_sc
200 integer iretcode,loc_nfun,orig_maxfun,n_try
201 double precision sc_dist,sc_dist_cutoff
202 double precision energy(0:n_ene),orig_e,cur_e
203 double precision evdw,escloc
204 double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
205 double precision var(maxvar)
207 double precision orig_theta(1:nres),orig_phi(1:nres),
208 + orig_alph(1:nres),orig_omeg(1:nres)
211 c Define what is meant by "neighbouring side-chain"
214 c Don't do glycine or ends
216 if (i.eq.10 .or. i.eq.21) return
218 c Freeze everything (later will relax only selected side-chains)
226 c Find the neighbours of the side-chain to move
227 c and save initial variables
232 c Don't do glycine (itype(j)==10)
233 if (itype(i).ne.10) then
234 sc_dist=dist(nres+i,nres+res_pick)
236 sc_dist=sc_dist_cutoff
238 if (sc_dist.lt.sc_dist_cutoff) then
239 nres_moved=nres_moved+1
249 e_sc=wsc*evdw+wscloc*escloc
250 cd call etotal(energy)
251 cd print *,'new ',(energy(k),k=0,n_ene)
256 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
257 c Move the selected residue (don't worry if it fails)
258 call gen_side(itype(res_pick),theta(res_pick+1),
259 + alph(res_pick),omeg(res_pick),fail)
261 c Minimize the side-chains starting from the new arrangement
262 call geom_to_var(nvar,var)
267 crc orig_theta(i)=theta(i)
268 crc orig_phi(i)=phi(i)
269 crc orig_alph(i)=alph(i)
270 crc orig_omeg(i)=omeg(i)
273 call minimize_sc1(e_sc,var,iretcode,loc_nfun)
275 cv write(*,'(2i3,2f12.5,2i3)')
276 cv & res_pick,nres_moved,orig_e,e_sc-cur_e,
277 cv & iretcode,loc_nfun
279 c$$$ if (iretcode.eq.8) then
280 c$$$ write(iout,*)'Coordinates just after code 8'
283 c$$$ call flush(iout)
285 c$$$ theta(i)=orig_theta(i)
286 c$$$ phi(i)=orig_phi(i)
287 c$$$ alph(i)=orig_alph(i)
288 c$$$ omeg(i)=orig_omeg(i)
290 c$$$ write(iout,*)'Coordinates just before code 8'
293 c$$$ call flush(iout)
298 call var_to_geom(nvar,var)
300 c If a lower energy was found, update the current structure...
301 if (e_sc.lt.cur_e) then
303 cv call etotal(energy)
306 cd e_sc1=wsc*evdw+wscloc*escloc
307 cd print *,' new',e_sc1,energy(0)
308 cv print *,'new ',energy(0)
309 cd call enerprint(energy(0))
312 if (mask_side(i).eq.1) then
318 c ...else revert to the previous one
321 if (mask_side(i).eq.1) then
330 n_steps=n_steps+n_try
332 c Reset the minimization mask_r to false
338 c-------------------------------------------------------------
340 subroutine sc_minimize(etot,iretcode,nfun)
341 c Minimizes side-chains only, leaving backbone frozen
345 implicit real*8 (a-h,o-z)
348 include 'COMMON.CHAIN'
349 include 'COMMON.FFIELD'
352 double precision etot
353 integer iretcode,nfun
357 double precision orig_w(n_ene),energy(0:n_ene)
358 double precision var(maxvar)
361 c Set non side-chain weights to zero (minimization is faster)
362 c NOTE: e(2) does not actually depend on the side-chain, only CA
389 c Prepare to freeze backbone
396 c Minimize the side-chains
398 call geom_to_var(nvar,var)
399 call minimize(etot,var,iretcode,nfun)
400 call var_to_geom(nvar,var)
403 c Put the original weights back and calculate the full energy
424 c-------------------------------------------------------------
425 subroutine minimize_sc1(etot,x,iretcode,nfun)
426 implicit real*8 (a-h,o-z)
428 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
429 include 'COMMON.IOUNITS'
432 include 'COMMON.MINIM'
435 double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
436 double precision energia(0:n_ene)
437 external func,gradient,fdum
438 external func_restr1,grad_restr1
439 logical not_done,change,reduce
440 common /przechowalnia/ v
442 call deflt(2,iv,liv,lv,v)
443 * 12 means fresh start, dont call deflt
445 * max num of fun calls
446 if (maxfun.eq.0) maxfun=500
448 * max num of iterations
449 if (maxmin.eq.0) maxmin=1000
453 * selects output unit
456 * 1 means to print out result
458 * 1 means to print out summary stats
460 * 1 means to print initial x and d
462 * min val for v(radfac) default is 0.1
464 * max val for v(radfac) default is 4.0
467 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
468 * the sumsl default is 0.1
470 * false conv if (act fnctn decrease) .lt. v(34)
471 * the sumsl default is 100*machep
473 * absolute convergence
474 if (tolf.eq.0.0D0) tolf=1.0D-4
476 * relative convergence
477 if (rtolf.eq.0.0D0) rtolf=1.0D-4
479 * controls initial step size
481 * large vals of d correspond to small components of step
489 call x2xx(x,xx,nvar_restr)
490 call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
491 & iv,liv,lv,v,idum,rdum,fdum)
494 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
502 ************************************************************************
503 subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
504 implicit real*8 (a-h,o-z)
506 include 'COMMON.DERIV'
507 include 'COMMON.IOUNITS'
509 include 'COMMON.FFIELD'
510 include 'COMMON.INTERACT'
511 include 'COMMON.TIME1'
513 double precision energia(0:n_ene),evdw,escloc
515 double precision ufparm,e1,e2
524 c Intercept NaNs in the coordinates, before calling etotal
530 if (x_sum.ne.x_sum) then
531 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
538 call var_to_geom_restr(n,x)
541 cd write (iout,*) 'ETOTAL called from FUNC'
544 f=wsc*evdw+wscloc*escloc
545 cd call etotal(energia(0))
546 cd f=wsc*energia(1)+wscloc*energia(12)
547 cd print *,f,evdw,escloc,energia(0)
549 C Sum up the components of the Cartesian gradient.
553 gradx(j,i,icg)=wsc*gvdwx(j,i)
559 c-------------------------------------------------------
560 subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
561 implicit real*8 (a-h,o-z)
563 include 'COMMON.CHAIN'
564 include 'COMMON.DERIV'
566 include 'COMMON.INTERACT'
567 include 'COMMON.FFIELD'
568 include 'COMMON.IOUNITS'
571 double precision urparm(1)
572 dimension x(maxvar),g(maxvar)
575 if (nf-nfl+1) 20,30,40
576 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
577 c write (iout,*) 'grad 20'
580 30 call var_to_geom_restr(n,x)
583 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
587 C Convert the Cartesian gradient into internal-coordinate gradient.
593 IF (mask_phi(i+2).eq.1) THEN
598 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
599 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
612 IF (mask_theta(i+2).eq.1) THEN
618 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
619 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
629 if (itype(i).ne.10) then
630 IF (mask_side(i).eq.1) THEN
634 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
643 if (itype(i).ne.10) then
644 IF (mask_side(i).eq.1) THEN
648 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
656 C Add the components corresponding to local energy terms.
663 if (mask_phi(i).eq.1) then
665 g(ig)=g(ig)+gloc(igall,icg)
671 if (mask_theta(i).eq.1) then
673 g(ig)=g(ig)+gloc(igall,icg)
679 if (itype(i).ne.10) then
681 if (mask_side(i).eq.1) then
683 g(ig)=g(ig)+gloc(igall,icg)
690 cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
694 C-----------------------------------------------------------------------------
695 subroutine egb1(evdw)
697 C This subroutine calculates the interaction energy of nonbonded side chains
698 C assuming the Gay-Berne potential of interaction.
700 implicit real*8 (a-h,o-z)
704 include 'COMMON.LOCAL'
705 include 'COMMON.CHAIN'
706 include 'COMMON.DERIV'
707 include 'COMMON.NAMES'
708 include 'COMMON.INTERACT'
709 include 'COMMON.IOUNITS'
710 include 'COMMON.CALC'
711 include 'COMMON.CONTROL'
714 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
717 c if (icall.eq.0) lprn=.true.
727 dxi=dc_norm(1,nres+i)
728 dyi=dc_norm(2,nres+i)
729 dzi=dc_norm(3,nres+i)
730 dsci_inv=dsc_inv(itypi)
732 C Calculate SC interaction energy.
735 do j=istart(i,iint),iend(i,iint)
736 IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
739 dscj_inv=dsc_inv(itypj)
740 sig0ij=sigma(itypi,itypj)
741 chi1=chi(itypi,itypj)
742 chi2=chi(itypj,itypi)
749 alf12=0.5D0*(alf1+alf2)
750 C For diagnostics only!!!
763 dxj=dc_norm(1,nres+j)
764 dyj=dc_norm(2,nres+j)
765 dzj=dc_norm(3,nres+j)
766 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
768 C Calculate angle-dependent terms of energy and contributions to their
772 sig=sig0ij*dsqrt(sigsq)
773 rij_shift=1.0D0/rij-sig+sig0ij
774 C I hate to put IF's in the loops, but here don't have another choice!!!!
775 if (rij_shift.le.0.0D0) then
777 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
778 cd & restyp(itypi),i,restyp(itypj),j,
779 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
783 c---------------------------------------------------------------
784 rij_shift=1.0D0/rij_shift
786 e1=fac*fac*aa(itypi,itypj)
787 e2=fac*bb(itypi,itypj)
788 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
789 eps2der=evdwij*eps3rt
790 eps3der=evdwij*eps2rt
791 evdwij=evdwij*eps2rt*eps3rt
794 sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
795 epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
796 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
797 cd & restyp(itypi),i,restyp(itypj),j,
798 cd & epsi,sigm,chi1,chi2,chip1,chip2,
799 cd & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
800 cd & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
804 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
807 C Calculate gradient components.
808 e1=e1*eps1*eps2rt**2*eps3rt**2
809 fac=-expon*(e1+evdwij)*rij_shift
812 C Calculate the radial part of the gradient
816 C Calculate angular part of the gradient.
823 C-----------------------------------------------------------------------------