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
158 ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
159 call chainbuild_extconf
163 c-------------------------------------------------------------
165 subroutine single_sc_move(res_pick,n_maxtry,e_drop,
166 + n_steps,n_fun,e_sc)
167 c Perturb one side-chain (res_pick) and minimize the
168 c neighbouring region, keeping all CA's and non-neighbouring
170 c Try until e_drop energy improvement is achieved, or n_maxtry
171 c attempts have been made
172 c At the start, e_sc should contain the side-chain-only energy(0)
173 c nsteps and nfun for this move are ADDED to n_steps and n_fun
177 implicit real*8 (a-h,o-z)
180 include 'COMMON.INTERACT'
181 include 'COMMON.CHAIN'
182 include 'COMMON.MINIM'
183 include 'COMMON.FFIELD'
184 include 'COMMON.IOUNITS'
187 double precision dist
191 integer res_pick,n_maxtry
192 double precision e_drop
194 c Input/Output arguments
195 integer n_steps,n_fun
196 double precision e_sc
202 integer iretcode,loc_nfun,orig_maxfun,n_try
203 double precision sc_dist,sc_dist_cutoff
204 double precision energy(0:n_ene),orig_e,cur_e
205 double precision evdw,escloc
206 double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
207 double precision var(maxvar)
209 double precision orig_theta(1:nres),orig_phi(1:nres),
210 + orig_alph(1:nres),orig_omeg(1:nres)
213 c Define what is meant by "neighbouring side-chain"
216 c Don't do glycine or ends
218 if (i.eq.10 .or. i.eq.ntyp1) return
220 c Freeze everything (later will relax only selected side-chains)
228 c Find the neighbours of the side-chain to move
229 c and save initial variables
234 c Don't do glycine (itype(j)==10)
235 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
236 sc_dist=dist(nres+i,nres+res_pick)
238 sc_dist=sc_dist_cutoff
240 if (sc_dist.lt.sc_dist_cutoff) then
241 nres_moved=nres_moved+1
248 call chainbuild_extconf
251 e_sc=wsc*evdw+wscloc*escloc
252 c write (iout,*) "sc_move: e_sc",e_sc
253 cd call etotal(energy)
254 cd print *,'new ',(energy(k),k=0,n_ene)
259 do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
260 c Move the selected residue (don't worry if it fails)
261 call gen_side(iabs(itype(res_pick)),theta(res_pick+1),
262 + alph(res_pick),omeg(res_pick),fail)
264 c Minimize the side-chains starting from the new arrangement
265 call geom_to_var(nvar,var)
270 crc orig_theta(i)=theta(i)
271 crc orig_phi(i)=phi(i)
272 crc orig_alph(i)=alph(i)
273 crc orig_omeg(i)=omeg(i)
276 call minimize_sc1(e_sc,var,iretcode,loc_nfun)
277 c write (iout,*) "n_try",n_try
278 c write (iout,*) "sc_move after minimze_sc1 e_sc",e_sc
279 cv write(*,'(2i3,2f12.5,2i3)')
280 cv & res_pick,nres_moved,orig_e,e_sc-cur_e,
281 cv & iretcode,loc_nfun
283 c$$$ if (iretcode.eq.8) then
284 c$$$ write(iout,*)'Coordinates just after code 8'
287 c$$$ call flush(iout)
289 c$$$ theta(i)=orig_theta(i)
290 c$$$ phi(i)=orig_phi(i)
291 c$$$ alph(i)=orig_alph(i)
292 c$$$ omeg(i)=orig_omeg(i)
294 c$$$ write(iout,*)'Coordinates just before code 8'
297 c$$$ call flush(iout)
302 call var_to_geom(nvar,var)
304 c If a lower energy was found, update the current structure...
305 if (e_sc.lt.cur_e) then
307 cv call etotal(energy)
310 cd e_sc1=wsc*evdw+wscloc*escloc
311 cd print *,' new',e_sc1,energy(0)
312 cv print *,'new ',energy(0)
313 cd call enerprint(energy(0))
316 if (mask_side(i).eq.1) then
322 c ...else revert to the previous one
325 if (mask_side(i).eq.1) then
334 n_steps=n_steps+n_try
336 c Reset the minimization mask_r to false
341 c-------------------------------------------------------------
342 subroutine minimize_sc1(etot,x,iretcode,nfun)
350 implicit real*8 (a-h,o-z)
353 c parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
354 parameter(max_sc_move=10)
355 parameter (liv=60,lv=(77+2*max_sc_move*(2*max_sc_move+17)/2))
357 include 'COMMON.IOUNITS'
360 include 'COMMON.MINIM'
362 double precision x(maxvar),d(maxvar),xx(maxvar)
363 double precision energia(0:n_ene)
366 common /zmienne/ nvar_restr
367 double precision grdmin
368 double precision funcgrad_restr1
369 external funcgrad_restr1
373 external func_restr1,grad_restr1
374 logical not_done,change,reduce
376 double precision v(1:lv)
377 common /przechowalnia/ v
381 coordtype='RIGIDBODY'
384 c jprint=print_min_stat
387 if (.not. allocated(scale)) allocate (scale(nvar))
389 c set scaling parameter for function and derivative values;
390 c use square root of median eigenvalue of typical Hessian
392 call x2xx(x,xx,nvar_restr)
403 c write (iout,*) "Calling lbfgs"
404 call lbfgs (nvar_restr,xx,etot,grdmin,funcgrad_restr1,optsave)
406 c write (iout,*) "After lbfgs"
409 call deflt(2,iv,liv,lv,v)
410 * 12 means fresh start, dont call deflt
412 * max num of fun calls
413 if (maxfun.eq.0) maxfun=500
415 * max num of iterations
416 if (maxmin.eq.0) maxmin=1000
420 * selects output unit
423 * 1 means to print out result
425 * 1 means to print out summary stats
427 * 1 means to print initial x and d
429 * min val for v(radfac) default is 0.1
431 * max val for v(radfac) default is 4.0
434 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
435 * the sumsl default is 0.1
437 * false conv if (act fnctn decrease) .lt. v(34)
438 * the sumsl default is 100*machep
440 * absolute convergence
441 if (tolf.eq.0.0D0) tolf=1.0D-4
443 * relative convergence
444 if (rtolf.eq.0.0D0) rtolf=1.0D-4
446 * controls initial step size
448 * large vals of d correspond to small components of step
456 call x2xx(x,xx,nvar_restr)
457 call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
458 & iv,liv,lv,v,idum,rdum,fdum)
461 c call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
470 double precision function funcgrad_restr1(x,g)
471 implicit real*8 (a-h,o-z)
473 include 'COMMON.DERIV'
474 include 'COMMON.IOUNITS'
476 include 'COMMON.FFIELD'
477 include 'COMMON.INTERACT'
478 include 'COMMON.TIME1'
479 include 'COMMON.CHAIN'
482 common /zmienne/ nvar_restr
483 double precision energia(0:n_ene),evdw,escloc
484 double precision ufparm,e1,e2
485 dimension x(maxvar),g(maxvar),gg(maxvar)
487 c Intercept NaNs in the coordinates, before calling etotal
493 if (x_sum.ne.x_sum) then
494 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
502 if (isnan(x(i))) then
506 write (iout,*) "NaN in coordinates"
512 c write (iout,*) "nvar_restr",nvar_restr
513 c write (iout,*) "x",(x(i),i=1,nvar_restr)
514 call var_to_geom_restr(nvar_restr,x)
516 call chainbuild_extconf
517 cd write (iout,*) 'ETOTAL called from FUNC'
520 f=wsc*evdw+wscloc*escloc
521 c write (iout,*) "evdw",evdw," escloc",escloc
528 c write (iout,*) "f",f
529 cd call etotal(energia(0))
530 cd f=wsc*energia(1)+wscloc*energia(12)
531 cd print *,f,evdw,escloc,energia(0)
533 C Sum up the components of the Cartesian gradient.
537 gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i)
542 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
544 call cart2intgrad(nvar,gg)
546 C Convert the Cartesian gradient into internal-coordinate gradient.
551 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
552 IF (mask_side(i).eq.1) THEN
555 c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
556 c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
561 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
562 IF (mask_side(i).eq.1) THEN
564 g(ig)=gg(ialph(i,1)+nside)
565 c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
566 c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
572 C Add the components corresponding to local energy terms.
579 if (mask_phi(i).eq.1) then
581 g(ig)=g(ig)+gloc(igall,icg)
587 if (mask_theta(i).eq.1) then
589 g(ig)=g(ig)+gloc(igall,icg)
595 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
597 if (mask_side(i).eq.1) then
599 g(ig)=g(ig)+gloc(igall,icg)
600 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
601 c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
608 cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
613 ************************************************************************
614 subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
615 implicit real*8 (a-h,o-z)
617 include 'COMMON.DERIV'
618 include 'COMMON.IOUNITS'
620 include 'COMMON.FFIELD'
621 include 'COMMON.INTERACT'
622 include 'COMMON.TIME1'
623 double precision energia(0:n_ene),evdw,escloc
624 double precision ufparm,e1,e2
633 c Intercept NaNs in the coordinates, before calling etotal
639 if (x_sum.ne.x_sum) then
640 write(iout,*)" *** func_restr1 : Found NaN in coordinates"
647 call var_to_geom_restr(n,x)
649 call chainbuild_extconf
650 cd write (iout,*) 'ETOTAL called from FUNC'
653 f=wsc*evdw+wscloc*escloc
654 c write (iout,*) "f",f
655 cd call etotal(energia(0))
656 cd f=wsc*energia(1)+wscloc*energia(12)
657 cd print *,f,evdw,escloc,energia(0)
659 C Sum up the components of the Cartesian gradient.
663 gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i)
669 c-------------------------------------------------------
670 subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
671 implicit real*8 (a-h,o-z)
673 include 'COMMON.CHAIN'
674 include 'COMMON.DERIV'
676 include 'COMMON.INTERACT'
677 include 'COMMON.FFIELD'
678 include 'COMMON.IOUNITS'
681 double precision urparm(1)
682 dimension x(maxvar),g(maxvar),gg(maxvar)
685 if (nf-nfl+1) 20,30,40
686 20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
687 c write (iout,*) 'grad 20'
690 30 call var_to_geom_restr(n,x)
691 call chainbuild_extconf
693 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
695 40 call cart2intgrad(nvar,gg)
697 C Convert the Cartesian gradient into internal-coordinate gradient.
703 IF (mask_phi(i+2).eq.1) THEN
711 IF (mask_theta(i+2).eq.1) THEN
718 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
719 IF (mask_side(i).eq.1) THEN
722 c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
723 c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
730 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
731 IF (mask_side(i).eq.1) THEN
733 g(ig)=gg(ialph(i,1)+nside)
734 c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
735 c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
741 C Add the components corresponding to local energy terms.
748 if (mask_phi(i).eq.1) then
750 g(ig)=g(ig)+gloc(igall,icg)
756 if (mask_theta(i).eq.1) then
758 g(ig)=g(ig)+gloc(igall,icg)
764 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
766 if (mask_side(i).eq.1) then
768 g(ig)=g(ig)+gloc(igall,icg)
769 c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
770 c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
777 cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
782 C-----------------------------------------------------------------------------
783 subroutine egb1(evdw)
785 C This subroutine calculates the interaction energy of nonbonded side chains
786 C assuming the Gay-Berne potential of interaction.
788 implicit real*8 (a-h,o-z)
792 include 'COMMON.LOCAL'
793 include 'COMMON.CHAIN'
794 include 'COMMON.DERIV'
795 include 'COMMON.NAMES'
796 include 'COMMON.INTERACT'
797 include 'COMMON.IOUNITS'
798 include 'COMMON.CALC'
799 include 'COMMON.CONTROL'
802 c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
805 c if (icall.eq.0) lprn=.true.
807 c do i=iatsc_s,iatsc_e
812 if (itypi.eq.ntyp1 .or. mask_side(i).eq.0) cycle
813 itypi1=iabs(itype(i+1))
818 if (xi.lt.0) xi=xi+boxxsize
820 if (yi.lt.0) yi=yi+boxysize
822 if (zi.lt.0) zi=zi+boxzsize
823 if ((zi.gt.bordlipbot)
824 &.and.(zi.lt.bordliptop)) then
825 C the energy transfer exist
826 if (zi.lt.buflipbot) then
827 C what fraction I am in
829 & ((positi-bordlipbot)/lipbufthick)
830 C lipbufthick is thickenes of lipid buffore
831 sslipi=sscalelip(fracinbuf)
832 ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
833 elseif (zi.gt.bufliptop) then
834 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
835 sslipi=sscalelip(fracinbuf)
836 ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
846 dxi=dc_norm(1,nres+i)
847 dyi=dc_norm(2,nres+i)
848 dzi=dc_norm(3,nres+i)
849 dsci_inv=dsc_inv(itypi)
851 C Calculate SC interaction energy.
853 c do iint=1,nint_gr(i)
854 c do j=istart(i,iint),iend(i,iint)
856 IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
859 if (itypj.eq.ntyp1) cycle
860 dscj_inv=dsc_inv(itypj)
861 sig0ij=sigma(itypi,itypj)
862 chi1=chi(itypi,itypj)
863 chi2=chi(itypj,itypi)
870 alf12=0.5D0*(alf1+alf2)
871 C For diagnostics only!!!
888 if (xj.lt.0) xj=xj+boxxsize
890 if (yj.lt.0) yj=yj+boxysize
892 if (zj.lt.0) zj=zj+boxzsize
893 if ((zj.gt.bordlipbot)
894 &.and.(zj.lt.bordliptop)) then
895 C the energy transfer exist
896 if (zj.lt.buflipbot) then
897 C what fraction I am in
899 & ((positi-bordlipbot)/lipbufthick)
900 C lipbufthick is thickenes of lipid buffore
901 sslipj=sscalelip(fracinbuf)
902 ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
903 elseif (zi.gt.bufliptop) then
904 fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
905 sslipj=sscalelip(fracinbuf)
906 ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
915 aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
916 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
917 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
918 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
920 dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
928 xj=xj_safe+xshift*boxxsize
929 yj=yj_safe+yshift*boxysize
930 zj=zj_safe+zshift*boxzsize
931 dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
932 if(dist_temp.lt.dist_init) then
942 if (subchap.eq.1) then
952 dxj=dc_norm(1,nres+j)
953 dyj=dc_norm(2,nres+j)
954 dzj=dc_norm(3,nres+j)
955 rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
957 sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
958 sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
960 C Calculate angle-dependent terms of energy and contributions to their
964 sig=sig0ij*dsqrt(sigsq)
965 rij_shift=1.0D0/rij-sig+sig0ij
966 C I hate to put IF's in the loops, but here don't have another choice!!!!
967 if (rij_shift.le.0.0D0) then
969 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
970 cd & restyp(itypi),i,restyp(itypj),j,
971 cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
975 c---------------------------------------------------------------
976 rij_shift=1.0D0/rij_shift
980 evdwij=eps1*eps2rt*eps3rt*(e1+e2)
981 eps2der=evdwij*eps3rt
982 eps3der=evdwij*eps2rt
983 evdwij=evdwij*eps2rt*eps3rt
986 sigm=dabs(aa/bb)**(1.0D0/6.0D0)
988 cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
989 cd & restyp(itypi),i,restyp(itypj),j,
990 cd & epsi,sigm,chi1,chi2,chip1,chip2,
991 cd & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
992 cd & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
996 if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
999 C Calculate gradient components.
1000 e1=e1*eps1*eps2rt**2*eps3rt**2
1001 fac=-expon*(e1+evdwij)*rij_shift
1004 fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
1005 C Calculate the radial part of the gradient
1009 gg_lipi(3)=ssgradlipi*evdwij
1010 gg_lipj(3)=ssgradlipj*evdwij
1011 C Calculate angular part of the gradient.
1018 C-----------------------------------------------------------------------------