double precision orig_w(n_ene)
double precision wtime
-
+ sideonly=.true.
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
wtor=orig_w(13)
wtor_d=orig_w(14)
wvdwpp=orig_w(15)
-
+ sideonly=.false.
+ mask_side=1
crc n_fun=n_fun+1
ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
return
nres_moved=0
do i=2,nres-1
c Don't do glycine (itype(j)==10)
- if (itype(i).ne.10) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
sc_dist=dist(nres+i,nres+res_pick)
else
sc_dist=sc_dist_cutoff
endif
enddo
- call chainbuild
+ call chainbuild_extconf
call egb1(evdw)
call esc(escloc)
e_sc=wsc*evdw+wscloc*escloc
+c write (iout,*) "sc_move: e_sc",e_sc
cd call etotal(energy)
cd print *,'new ',(energy(k),k=0,n_ene)
orig_e=e_sc
crc enddo
call minimize_sc1(e_sc,var,iretcode,loc_nfun)
-
+c write (iout,*) "n_try",n_try
+c write (iout,*) "sc_move after minimze_sc1 e_sc",e_sc
cv write(*,'(2i3,2f12.5,2i3)')
cv & res_pick,nres_moved,orig_e,e_sc-cur_e,
cv & iretcode,loc_nfun
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)
+#ifdef LBFGS_SC
+ use minima
+ use inform
+ use output
+ use iounit
+ use scales
+#endif
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+#ifndef LBFGS_SC
+c parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+ parameter(max_sc_move=10)
+ parameter (liv=60,lv=(77+2*max_sc_move*(2*max_sc_move+17)/2))
+#endif
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 x(maxvar),d(maxvar),xx(maxvar)
double precision energia(0:n_ene)
+#ifdef LBFGS_SC
+ integer nvar_restr
+ common /zmienne/ nvar_restr
+ double precision grdmin
+ double precision funcgrad_restr1
+ external funcgrad_restr1
+ external optsave
+#else
external func,gradient,fdum
external func_restr1,grad_restr1
logical not_done,change,reduce
+ dimension iv(liv)
+ double precision v(1:lv)
common /przechowalnia/ v
-
+#endif
+#ifdef LBFGS_SC
+ maxiter=7
+ coordtype='RIGIDBODY'
+ grdmin=tolf
+ jout=iout
+c jprint=print_min_stat
+ jprint=0
+ iwrite=0
+ if (.not. allocated(scale)) allocate (scale(nvar))
+c
+c set scaling parameter for function and derivative values;
+c use square root of median eigenvalue of typical Hessian
+c
+ call x2xx(x,xx,nvar_restr)
+ set_scale = .true.
+c nvar = 0
+ do i = 1, nvar_restr
+c if (use(i)) then
+c do j = 1, 3
+c nvar = nvar + 1
+ scale(i) = 12.0d0
+c end do
+c end if
+ end do
+c write (iout,*) "Calling lbfgs"
+ call lbfgs (nvar_restr,xx,etot,grdmin,funcgrad_restr1,optsave)
+ deallocate(scale)
+c write (iout,*) "After lbfgs"
+ call xx2x(x,xx)
+#else
call deflt(2,iv,liv,lv,v)
* 12 means fresh start, dont call deflt
iv(1)=12
* controls output
iv(19)=2
* selects output unit
-c iv(21)=iout
iv(21)=0
+c iv(21)=0
* 1 means to print out result
iv(22)=0
* 1 means to print out summary stats
& 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)
+c call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
ENDIF
etot=v(10)
iretcode=iv(1)
nfun=iv(6)
-
+#endif
return
end
+#ifdef LBFGS_SC
+ double precision function funcgrad_restr1(x,g)
+ 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'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ integer nvar_restr
+ common /zmienne/ nvar_restr
+ double precision energia(0:n_ene),evdw,escloc
+ double precision ufparm,e1,e2
+ dimension x(maxvar),g(maxvar),gg(maxvar)
+#ifdef OSF
+c Intercept NaNs in the coordinates, before calling etotal
+ x_sum=0.D0
+ do i=1,nvar_restr
+ 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
+#else
+ FOUND_NAN=.false.
+ do i=1,nvar_restr
+ if (isnan(x(i))) then
+ FOUND_NAN=.true.
+ f=1.0D+73
+ funcgrad_restr1=f
+ write (iout,*) "NaN in coordinates"
+ return
+ endif
+ enddo
+#endif
+
+c write (iout,*) "nvar_restr",nvar_restr
+c write (iout,*) "x",(x(i),i=1,nvar_restr)
+ call var_to_geom_restr(nvar_restr,x)
+ call zerograd
+ call chainbuild_extconf
+cd write (iout,*) 'ETOTAL called from FUNC'
+ call egb1(evdw)
+ call esc(escloc)
+ f=wsc*evdw+wscloc*escloc
+c write (iout,*) "evdw",evdw," escloc",escloc
+ if (isnan(f)) then
+ f=1.0d20
+ funcgrad_restr1=f
+ return
+ endif
+ funcgrad_restr1=f
+c write (iout,*) "f",f
+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)+wscloc*gsclocx(j,i)
+ enddo
+ enddo
+
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+ call cart2intgrad(nvar,gg)
+C
+C Convert the Cartesian gradient into internal-coordinate gradient.
+C
+
+ ig=0
+ do i=2,nres-1
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ IF (mask_side(i).eq.1) THEN
+ ig=ig+1
+ g(ig)=gg(ialph(i,1))
+c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
+c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
+ ENDIF
+ endif
+ enddo
+ do i=2,nres-1
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ IF (mask_side(i).eq.1) THEN
+ ig=ig+1
+ g(ig)=gg(ialph(i,1)+nside)
+c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
+c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
+ 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 .and. itype(i).ne.ntyp1) then
+ igall=igall+1
+ if (mask_side(i).eq.1) then
+ ig=ig+1
+ g(ig)=g(ig)+gloc(igall,icg)
+c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
+c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
+ 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
+#else
************************************************************************
subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
implicit real*8 (a-h,o-z)
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)
call var_to_geom_restr(n,x)
call zerograd
- call chainbuild
+ call chainbuild_extconf
cd write (iout,*) 'ETOTAL called from FUNC'
call egb1(evdw)
call esc(escloc)
f=wsc*evdw+wscloc*escloc
+c write (iout,*) "f",f
cd call etotal(energia(0))
cd f=wsc*energia(1)+wscloc*energia(12)
cd print *,f,evdw,escloc,energia(0)
C
do i=1,nct
do j=1,3
- gradx(j,i,icg)=wsc*gvdwx(j,i)
+ gradx(j,i,icg)=wsc*gvdwx(j,i)+wscloc*gsclocx(j,i)
enddo
enddo
external ufparm
integer uiparm(1)
double precision urparm(1)
- dimension x(maxvar),g(maxvar)
+ dimension x(maxvar),g(maxvar),gg(maxvar)
icg=mod(nf,2)+1
if (nf-nfl+1) 20,30,40
if (nf.eq.0) return
goto 40
30 call var_to_geom_restr(n,x)
- call chainbuild
+ call chainbuild_extconf
C
C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
C
- 40 call cartder
+ 40 call cart2intgrad(nvar,gg)
C
C Convert the Cartesian gradient into internal-coordinate gradient.
C
ig=0
- ind=nres-2
+ 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
+ IF (mask_phi(i+2).eq.1) THEN
ig=ig+1
- g(ig)=gphii
- ELSE
- ind=ind+nres-1-i
+ g(ig)=gg(i-1)
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
+ g(ig)=gg(nphi+i)
ENDIF
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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
+ g(ig)=gg(ialph(i,1))
+c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)
+c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1))
ENDIF
endif
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) 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
+ g(ig)=gg(ialph(i,1)+nside)
+c write (iout,*) "i",i," ig",ig," ialph",ialph(i,1)+nside
+c write (iout,*) "g",g(ig)," gg",gg(ialph(i,1)+nside)
ENDIF
endif
enddo
do ij=1,2
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
igall=igall+1
if (mask_side(i).eq.1) then
ig=ig+1
g(ig)=g(ig)+gloc(igall,icg)
+c write (iout,*) "ij",ij," i",i," ig",ig," igall",igall
+c write (iout,*) "gloc",gloc(igall,icg)," g",g(ig)
endif
endif
enddo
cd enddo
return
end
+#endif
C-----------------------------------------------------------------------------
subroutine egb1(evdw)
C
lprn=.false.
c if (icall.eq.0) lprn=.true.
ind=0
- do i=iatsc_s,iatsc_e
+c do i=iatsc_s,iatsc_e
+ do i=nnt,nct
itypi=iabs(itype(i))
- if (itypi.eq.ntyp1) cycle
+ if (itypi.eq.ntyp1 .or. mask_side(i).eq.0) cycle
itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
C
C Calculate SC interaction energy.
C
- do iint=1,nint_gr(i)
- do j=istart(i,iint),iend(i,iint)
+c do iint=1,nint_gr(i)
+c do j=istart(i,iint),iend(i,iint)
+ do j=i+1,nct
IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
ind=ind+1
itypj=iabs(itype(j))
call sc_grad
ENDIF
enddo ! j
- enddo ! iint
+c enddo ! iint
enddo ! i
end
C-----------------------------------------------------------------------------
projects / unres.git / blobdiff