data/MPI_data.F90
data/REMD_data.F90
data/comm_local.F90
+ inform.F90
+ keys.F90
+ math2.F90
+ iounit.F90
+ linmin.F90
+ minima.F90
+ output.F90
+ scales.F90
+ search.F90
+ optsave_dum.F90
prng_32.F90
math.F90
random.F90
geometry.F90
- io_base.F90
+ io_base.F90
energy.F90
check_bond.F90
control.F90
set (CMAKE_Fortran_FLAGS_RELEASE " ")
set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -traceback")
# set(FFLAGS0 "-fpp -c -CB -g -ip " )
- set(FFLAGS0 "-O3 -g -ip -fpp -mcmodel=large" )
+ set(FFLAGS0 "-CB -g -ip -fpp -mcmodel=large" )
# set(FFLAGS0 "-O0 -CB -CA -g" )
set(FFLAGS1 "-fpp -c -O " )
set(FFLAGS2 "-fpp -c -g -CA -CB ")
if(UNRES_NEWGRAD)
# set(CPPFLAGS "${CPPFLAGS} -DFIVEDIAG ")
-set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG -DSC_END")
+set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG -DSC_END -DLBFGS")
endif()
# set(CPPFLAGS "${CPPFLAGS} -DCARGRAD -DFIVEDIAG ")
!-----------------------------------------------------------------------------
! Maximum number of moves (n1-n8)
integer,parameter :: mxmv=18
+
!-----------------------------------------------------------------------------
!
!
! minim_jlee.F
!-----------------------------------------------------------------------------
subroutine minim_jlee
-
+#ifdef LBFGS
+ use minima
+ use inform
+ use output
+ use iounit
+ use scales
+ use energy, only: funcgrad
+#endif
use minim_data
use MPI_data
use energy_data
use energy, only:fdum
use control, only:init_int_table
use minimm, only:sumsl,deflt
+#ifdef LBFGS
+ use minimm,only: lbfgs
+#endif
! controls minimization and sorting routines
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
include 'mpif.h'
integer,parameter :: liv=60
integer :: lv
+#ifdef LBFGS
+ double precision grdmin
+ real(kind=8),dimension(6*nres)::g
+ external optsave
+#endif
! external func,gradient!,fdum !use minim & energy
! real(kind=4) :: ran1,ran2,ran3
! include 'COMMON.SETUP'
logical :: fail !check_var,
integer :: iloop(2)
!el common /przechowalnia/ v
- integer :: i,j,ierr,n,nfun,nft_sc,nf,ierror,ierrcode
+ integer :: i,j,ierr,n,nft_sc,nf,ierror,ierrcode
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: rad,eee,etot !,fdum
!el from subroutine parmread
! Define the constants of the disulfide bridge
! receive # of start
! print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
! & ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
+#ifndef LBFGS
if (.not. allocated(v)) allocate(v(1:lv))
+#endif
+#ifdef LBFGS
+ maxiter=maxmin
+ coordtype='RIGIDBODY'
+ grdmin=tolf
+ jout=iout
+ jprint=print_min_stat
+ 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
+ set_scale = .true.
+!c nvar = 0
+ do i = 1, nvar
+!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
+#endif
+
nhpb0=nhpb
10 continue
time0s=MPI_WTIME()
nfun=nfun+1
write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
else
+! v(10)=1.0d20
+! iv(1)=-1
+#ifdef LBFGS
+ etot=1.0d20
+ nfun=-1
+#else
v(10)=1.0d20
iv(1)=-1
+#endif
goto 201
endif
endif
endif
if (check_var(var,info)) then
+! v(10)=1.0d21
+! iv(1)=6
+#ifdef LBFGS
+ etot=1.0d21
+#else
v(10)=1.0d21
iv(1)=6
+#endif
+
goto 201
endif
do i=1,nvar
garbage(i)=var(i)
enddo
+#ifdef LBFGS
+ eee=funcgrad(var,g)
+ nfun=nfun+1
+ if(eee.ge.1.0d20) then
+!c print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
+!c print *,' energy before SUMSL =',eee
+!c print *,' aborting local minimization'
+ go to 201
+ endif
+
+ call lbfgs (nvar,var,etot,grdmin,funcgrad,optsave)
+ deallocate(scale)
+#else
call deflt(2,iv,liv,lv,v)
! 12 means fresh start, dont call deflt
! find which conformation was returned from sumsl
nfun=nfun+iv(7)
+#endif
! print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
! & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
! if (iv(1).ne.4 .or. nf.le.1) then
201 continue
indx(1)=n
! return code: 6-gradient 9-number of ftn evaluation, etc
+! indx(2)=iv(1)
+#ifdef LBFGS
+ indx(2)=nfun
+#else
indx(2)=iv(1)
+#endif
! total # of ftn evaluations (for iwf=0, it includes all minimizations).
indx(3)=nfun
indx(4)=info(2)
! send back energies
! al & cc
! calculate contact order
+#ifdef LBFGS
#ifdef CO_BIAS
- call contact(.false.,ncont,icont_,co)
+ call contact(.false.,nnt,nct,ncont,icont,co)
+ erg(1)=etot-1.0d2*co
+#else
+ erg(1)=etot
+#endif
+#else
+#ifdef CO_BIAS
+ call contact(.false.,nnt,nct,ncont,icont,co)
erg(1)=v(10)-1.0d2*co
#else
erg(1)=v(10)
#endif
+#endif
j=1
call mpi_send(erg,j,mpi_double_precision,king,idreal,&
CG_COMM,ierr)
real(kind=8),dimension(0:n_ene) :: energia,energia_ave
!!! local variables -el
- integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,&
+ integer :: i,ii,kkk,it,j,nacc,ijunk,indmin,indmax,&
ISWEEP,Kwita,iretcode,indeold,iene,noverlap,&
irep,nstart_grow,inde
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,&
deix,dent,sold,scur,runtime
!
real(kind=8),dimension(0:n_ene) :: energia,energia_ave
!
!!! local variables - el
- integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,&
+ integer :: i,j,it,ii,iproc,nacc,ISWEEP,indmax,indmin,ijunk,&
Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,&
frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms
real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres
!-----------------------------------------------------------------------------
! common /przechowalnia/ subroutine: setup_fricmat
+#ifndef LBFGS
real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres)
+#endif
!-----------------------------------------------------------------------------
!
!
! include 'COMMON.NAMES'
! include 'COMMON.TIME1'
! include 'COMMON.HAIRPIN'
- real(kind=8),dimension(3) :: L,vcm
+ real(kind=8),dimension(3) :: L,vcm,boxx
#ifdef VOUT
real(kind=8),dimension(6*nres) :: v_work,v_transf !(maxres6) maxres6=6*maxres
#endif
!el common /gucio/ cm
integer :: i,j,nharp
integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
+
! logical :: ovrtim
real(kind=8) :: tt0,scalfac
integer :: nres2,itime
nres2=2*nres
print *, "ENTER MD"
+ boxx(1)=boxxsize
+ boxx(2)=boxysize
+ boxx(3)=boxzsize
+
!
#ifdef MPI
print *,"MY tmpdir",tmpdir,ilen(tmpdir)
endif
if (reset_vel .and. tbf .and. lang.eq.0 &
.and. mod(itime,count_reset_vel).eq.0) then
+ !WARP WATER
+ do i=1,nres
+ if (molnum(i).eq.5) then
+ call to_box(c(1,i),c(2,i),c(3,i))
+ do j=1,3
+ if (c(j,i).le.0) c(j,i)=c(j,i)+boxx(j)
+ enddo
+ endif
+! write(iout,*) "COORD",c(1,i),c(2,i),c(3,i)
+ enddo
call random_vel
+
write(iout,'(a,f20.2)') &
"Velocities reset to random values, time",totT
do i=0,2*nres
character(len=50) :: tytul
logical :: file_exist
!el common /gucio/ cm
- integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,ierr,mnum,itime
+ integer :: i,j,ipos,iq,iw,nft_sc,iretcode,ierr,mnum,itime
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: etot,tt0
logical :: fail
if (me.eq.king.or..not.out1file) write(iout,*)&
'Minimizing initial PDB structure: Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)&
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
+#endif
else
call geom_to_var(nvar,varia)
if(me.eq.king.or..not.out1file) write (iout,*)&
call var_to_geom(nvar,varia)
#ifdef LBFGS
if (me.eq.king.or..not.out1file)&
- write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
if(me.eq.king.or..not.out1file)&
- write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ write(iout,*) 'LBFGS return code is ',statusbf,' eval ',nfun
#else
if (me.eq.king.or..not.out1file)&
write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
!#define DEBUG
#ifdef FIVEDIAG
real(kind=8) ,allocatable, dimension(:) :: xsolv,DML,rs
- real(kind=8) :: sumx,Ek2,Ek3,aux
+ real(kind=8) :: sumx,Ek2,Ek3,aux,masinv
#ifdef DEBUG
real(kind=8) ,allocatable, dimension(:) :: rsold
real (kind=8),allocatable,dimension(:,:) :: matold,inertia
#endif
do ichain=1,nchain
ind=0
- if(.not.allocated(ghalf)) print *,"COCO"
- if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
- ghalf=0.0d0
+! if(.not.allocated(ghalf)) print *,"COCO"
+! if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2))
+! ghalf=0.0d0
n=dimen_chain(ichain)
innt=iposd_chain(ichain)
+ if (molnum(innt).eq.5) go to 137
+ if(.not.allocated(ghalf)) print *,"COCO"
+ if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2))
+ ghalf=0.0d0
inct=innt+n-1
#ifdef DEBUG
write (iout,*) "Chain",ichain," n",n," start",innt
enddo
enddo
#endif
+137 continue
+ write(iout,*) "HERE,",n
xv=0.0d0
ii=0
do i=1,n
do k=1,3
ii=ii+1
+ if (molnum(innt).eq.5) geigen(i)=1.0/msc(itype(innt+i-1,5),5)
+
sigv=dsqrt((Rb*t_bath)/geigen(i))
lowb=-5*sigv
highb=5*sigv
d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
-!c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
-!c & " d_t_work_new",d_t_work_new(ii)
+! write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
+! " d_t_work_new",d_t_work_new(ii)
enddo
enddo
+ if (molnum(innt).eq.5) then
+
+ do k=1,3
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ masinv=1.0/msc(itype(innt+i-1,5),5)
+ d_t_work(ind)=d_t_work(ind)&
+ +masinv*d_t_work_new((i-1)*3+k)
+ enddo
+ enddo
+
+ else
do k=1,3
do i=1,n
ind=(i-1)*3+k
d_t_work(ind)=d_t_work(ind)&
+Gvec(i,j)*d_t_work_new((j-1)*3+k)
enddo
-!c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
-!c call flush(iout)
enddo
enddo
+ endif
#ifdef DEBUG
aux=0.0d0
do k=1,3
allocate(DM(nres2),DU1(nres2),DU2(nres2))
allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
!#ifdef DEBUG
- allocate(Gvec(nres2,nres2))
+ allocate(Gvec(1300,1300))!maximum number of protein in one chain
!#endif
#else
write (iout,*) "Before A Allocation",nfgtasks-1
ncont=0
kkk=3
- do i=nnt+kkk,nct
+ do i=nnt_molec(1)+kkk,nct_molec(1)
iti=iabs(itype(i,molnum(i)))
if (molnum(i).lt.3) then
inum=i+nres
! print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
if (dist(inum,jnum).lt.rcomp) then
ncont=ncont+1
+ if (ncont.gt.nres*100) ncont=nres*100
icont(1,ncont)=i
icont(2,ncont)=j
endif
ncont=0
kkk=0
! print *,'nnt=',nnt,' nct=',nct
- do i=nnt,nct-3
+ do i=nnt_molec(1),nct_molec(1)-3
do k=1,3
c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
enddo
ees=0.0
evdw=0.0
! print *, "nntt,nct",nnt,nct-2
- do 1 i=nnt,nct-2
+ do 1 i=nnt_molec(1),nct_molec(1)-2
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) goto 1
xi=c(1,i)
yi=c(2,i)
jgrad_start(i)=i+1
jgrad_end(i)=nres
enddo
+! THIS SHOULD BE FIXED
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+ call int_bounds(nres_molec(5),icatb_start,icatb_end)
+ icatb_start=icatb_start+itmp
+ icatb_end=icatb_end+itmp
+
+
+
if (lprint) then
write (*,*) 'Processor:',fg_rank,' CG group',kolor,&
' absolute rank',myrank,&
iphi_nucl_end,iphid_nucl_start,iphid_nucl_end,&
ibondp_nucl_start,ibondp_nucl_end,ithet_nucl_start,ithet_nucl_end,&
loc_start_nucl,loc_end_nucl
+ integer :: icatb_start,icatb_end
integer,dimension(:),allocatable :: ibond_displ,ibond_count,&
ithet_displ,ithet_count,iphi_displ,iphi_count,iphi1_displ,&
iphi1_count,ivec_displ,ivec_count,iset_displ,iset_count,&
newcontlistcatscangfi,newcontlistcatscangfj,&
newcontlistcatscangfk,&
newcontlistcatscangti,newcontlistcatscangtj,&
- newcontlistcatscangtk,newcontlistcatscangtl
+ newcontlistcatscangtk,newcontlistcatscangtl,&
+ newcontlistcatcatnormi,newcontlistcatcatnormj
g_ilist_catptran,g_ilist_catscang,g_ilist_catscangf,g_ilist_catscangt,&
g_listcatscang_start,g_listcatscang_end,&
g_listcatscangf_start,g_listcatscangf_end,&
- g_listcatscangt_start,g_listcatscangt_end
+ g_listcatscangt_start,g_listcatscangt_end,&
+ g_listcatcatnorm_start,g_listcatcatnorm_end,g_ilist_catcatnorm
! MARTINI FORCE FIELD
! buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
!-------------------------------------------------------------------------
real(kind=8),dimension(3,70000) :: ea
+#ifdef LBFGS
+ character*9 statusbf
+ integer niter,nfun,ncalls
+#endif
! real(kind=8) :: buftubebot, buftubetop,bordtubebot,bordtubetop,tubebufthick
end module energy_data
integer ishield_listbuf(-1:nres_molec(1)), &
shield_listbuf(maxcontsshi,-1:nres_molec(1)),k,j,i,iii,impishi,mojint,jjj
+ integer :: imatupdate2
! print *,"I START ENERGY"
imatupdate=100
+ imatupdate2=100
! if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! real(kind=8), dimension(:),allocatable:: fac_shieldbuf
! real(kind=8), dimension(:,:,:),allocatable:: &
time_Bcastw=time_Bcastw+MPI_Wtime()-time00
! call chainbuild_cart
endif
- print *,"itime_mat",itime_mat,imatupdate
+! print *,"itime_mat",itime_mat,imatupdate
if (nfgtasks.gt.1) then
call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
endif
if (mod(itime_mat,imatupdate).eq.0) then
! print *,'Processor',myrank,' calling etotal ipot=',ipot
call make_cat_pep_list
+! call make_cat_cat_list
endif
endif
endif
+ if (nres_molec(5).gt.0) then
+ if (mod(itime_mat,imatupdate2).eq.0) then
+! print *, "before cat cat"
+! print *,'Processor',myrank,' calling etotal ipot=',ipot
+! call make_cat_pep_list
+ call make_cat_cat_list
+ endif
+ endif
! write (iout,*) "after make_pp_inter_list"
! print *,'Processor',myrank,' calling etotal ipot=',ipot
Eafmforce=0.0d0
endif
endif
- print *,"before tubemode",tubemode
+! print *,"before tubemode",tubemode
if (tubemode.eq.1) then
call calctube(etube)
else if (tubemode.eq.2) then
else
etube=0.0d0
endif
- print *, "TU JEST PRZED EHPB"
+! print *, "TU JEST PRZED EHPB"
call edis(ehpb)
!--------------------------------------------------------
else
ecation_protang=0.0d0
endif
- if (nfgtasks.gt.1) then
- if (fg_rank.eq.0) then
+! if (nfgtasks.gt.1) then
+! if (fg_rank.eq.0) then
if (nres_molec(5).gt.1) call ecatcat(ecationcation)
- endif
- else
- if (nres_molec(5).gt.1) call ecatcat(ecationcation)
- endif
+! endif
+! else
+! if (nres_molec(5).gt.1) call ecatcat(ecationcation)
+! endif
if (oldion.gt.0) then
if (g_ilist_catpnorm.gt.0) call ecat_prot(ecation_prot)
else
! include 'COMMON.VAR'
! include 'COMMON.CONTACTS'
use comm_srutu
+!#ifdef LBFGS
+! use minimm, only: funcgrad
+!#endif
!el integer :: icall
!el common /srutu/ icall
+! real(kind=8) :: funcgrad
real(kind=8),dimension(6) :: ggg
real(kind=8),dimension(3) :: cc,xx,ddc,ddx
real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
real(kind=8) :: urparm(1)
!EL external fdum
integer :: nf,i,j,k
- real(kind=8) :: aincr,etot,etot1
+ real(kind=8) :: aincr,etot,etot1,ff
icg=1
nf=0
nfl=0
call geom_to_var(nvar,x)
call etotal(energia)
etot=energia(0)
+#ifdef LBFGS
+ ff=funcgrad(x,g)
+#else
!el call enerprint(energia)
call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+#endif
icall =1
do i=1,nres
write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
! include 'COMMON.VAR'
! include 'COMMON.GEO'
use comm_srutu
+!#ifdef LBFGS
+! use minimm, only : funcgrad
+!#endif
!el integer :: icall
!el common /srutu/ icall
+! real(kind=8) :: funcgrad
real(kind=8),dimension(6*nres) :: x,gana,gg !(maxvar) (maxvar=6*maxres)
integer :: uiparm(1)
real(kind=8) :: urparm(1)
character(len=6) :: key
!EL external fdum
integer :: i,ii,nf
- real(kind=8) :: xi,aincr,etot,etot1,etot2
+ real(kind=8) :: xi,aincr,etot,etot1,etot2,ff
call zerograd
aincr=1.0D-7
print '(a)','Calling CHECK_INT.'
#endif
nf=1
nfl=3
+#ifdef LBFGS
+ ff=funcgrad(x,gana)
+#else
+
!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
+#endif
icall=1
do i=1,nvar
xi=x(i)
!-----------------------------------------------------------------------------
! gradient_p.F
!-----------------------------------------------------------------------------
+#ifndef LBFGS
subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
use io_base, only:intout,briefout
!d write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
return
end subroutine gradient
+#endif
!-----------------------------------------------------------------------------
subroutine func(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
allocate(newcontlistcatscnormj(300*nres))
allocate(newcontlistcatpnormi(300*nres))
allocate(newcontlistcatpnormj(300*nres))
+ allocate(newcontlistcatcatnormi(900*nres))
+ allocate(newcontlistcatcatnormj(900*nres))
allocate(newcontlistcatscangi(300*nres))
allocate(newcontlistcatscangj(300*nres))
#endif
subroutine ecatcat(ecationcation)
use MD_data, only: t_bath
- integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff
+ integer :: i,j,itmp,xshift,yshift,zshift,subchap,k,itypi,itypj,irdiff,&
+ ii
real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
r7,r4,ecationcation,k0,rcal,aa,bb,sslipi,ssgradlipi,sslipj,ssgradlipj
real(kind=8) :: xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
! k0 = 332.0*(2.0*2.0)/80.0
itmp=0
- do i=1,4
- itmp=itmp+nres_molec(i)
- enddo
-! write(iout,*) "itmp",itmp
- do i=itmp+1,itmp+nres_molec(5)-1
-
+! do i=1,4
+! itmp=itmp+nres_molec(i)
+! enddo
+! write(iout,*) "itmp",g_listcatcatnorm_start, g_listcatcatnorm_end
+! do i=itmp+1,itmp+nres_molec(5)-1
+ do ii=g_listcatcatnorm_start, g_listcatcatnorm_end
+ i=newcontlistcatcatnormi(ii)
+ j=newcontlistcatcatnormj(ii)
+
xi=c(1,i)
yi=c(2,i)
zi=c(3,i)
itypi=itype(i,5)
call to_box(xi,yi,zi)
call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
- do j=i+1,itmp+nres_molec(5)
+! do j=i+1,itmp+nres_molec(5)
itypj=itype(j,5)
! print *,i,j,itypi,itypj
k0 = 332.0*(ichargecat(itypi)*ichargecat(itypj))/80.0
rcal =xj**2+yj**2+zj**2
ract=sqrt(rcal)
if ((itypi.gt.1).or.(itypj.gt.1)) then
-
+ if (sss2min2.eq.0.0d0) cycle
+ sss2min2=sscale2(ract,12.0d0,1.0d0)
+ sss2mingrad2=sscagrad2(ract,12.0d0,1.0d0)
! rcat0=3.472
! epscalc=0.05
! r06 = rcat0**6
enddo
do k=1,3
gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k)
- gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
- gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
+ gradcatcat(k,i)=gradcatcat(k,i)-(gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2)
+ gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+(Evan1cat+Evan2cat+Eeleccat)*sss2mingrad2
enddo
if (energy_dec) write (iout,*) "ecatcat",i,j,Evan1cat,Evan2cat,Eeleccat,&
r012,rcal**6,ichargecat(itypi)*ichargecat(itypj)
! write(iout,*) "ecatcat",i,j, ecationcation,xj,yj,zj
- ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
+ ecationcation=ecationcation+(Evan1cat+Evan2cat+Eeleccat)*sss2min2
else !this is water part and other non standard molecules
sss2min2=sscale2(ract,10.0d0,1.0d0)! cutoff for water interaction is 15A
gradcatcat(k,i)=gradcatcat(k,i)-gg(k)*sss2min2-sss2mingrad2*ewater*r(k)/ract
gradcatcat(k,j)=gradcatcat(k,j)+gg(k)*sss2min2+sss2mingrad2*ewater*r(k)/ract
enddo
- if (energy_dec) write(iout,'(2f10.7,f15.7,2i5)') rdiff,ract,ecationcation,i,j
+ if (energy_dec) write(iout,'(2f8.2,f10.2,2i5)') rdiff,ract,ecationcation,i,j
endif ! end water
enddo
- enddo
+! enddo
return
end subroutine ecatcat
!---------------------------------------------------------------------------
real(kind=8) :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
real(kind=8) :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
real(kind=8) :: xja,yja,zja
- integer:: contlistcatpnormi(250*nres),contlistcatpnormj(250*nres)
+ integer:: contlistcatpnormi(300*nres),contlistcatpnormj(300*nres)
integer:: contlistcatscnormi(250*nres),contlistcatscnormj(250*nres)
integer:: contlistcatptrani(250*nres),contlistcatptranj(250*nres)
integer:: contlistcatsctrani(250*nres),contlistcatsctranj(250*nres)
return
end subroutine make_cat_pep_list
+ subroutine make_cat_cat_list
+ include 'mpif.h'
+ real(kind=8) :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
+ real(kind=8) :: xmedj,ymedj,zmedj,sslipi,ssgradlipi,faclipij2,sslipj,ssgradlipj
+ real(kind=8) :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
+ real(kind=8) :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
+ real(kind=8) :: xja,yja,zja
+ integer,dimension(:),allocatable:: contlistcatpnormi,contlistcatpnormj
+! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
+ integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_catscnorm,&
+ ilist_catsctran,ilist_catpnorm,ilist_catptran,itmp,ilist_catscang,&
+ ilist_catscangf,ilist_catscangt,k
+ integer displ(0:nprocs),i_ilist_catscnorm(0:nprocs),ierr,&
+ i_ilist_catpnorm(0:nprocs),i_ilist_catsctran(0:nprocs),&
+ i_ilist_catptran(0:nprocs),i_ilist_catscang(0:nprocs),&
+ i_ilist_catscangf(0:nprocs),i_ilist_catscangt(0:nprocs)
+ write(iout,*),"START make_catcat"
+ ilist_catpnorm=0
+ ilist_catscnorm=0
+ ilist_catptran=0
+ ilist_catsctran=0
+ ilist_catscang=0
+
+ if (.not.allocated(contlistcatpnormi)) then
+ allocate(contlistcatpnormi(900*nres))
+ allocate(contlistcatpnormj(900*nres))
+ endif
+ r_buff_list=3.0
+ itmp=0
+ do i=1,4
+ itmp=itmp+nres_molec(i)
+ enddo
+! go to 17
+! do i=1,nres_molec(1)-1 ! loop over all peptide groups needs parralelization
+ do i=icatb_start,icatb_end
+ xi=c(1,i)
+ yi=c(2,i)
+ zi=c(3,i)
+ call to_box(xi,yi,zi)
+ dxi=dc_norm(1,i)
+ dyi=dc_norm(2,i)
+ dzi=dc_norm(3,i)
+! dsci_inv=vbld_inv(i+nres)
+ do j=i+1,itmp+nres_molec(5)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+! call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+! faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+ xja=boxshift(xj-xi,boxxsize)
+ yja=boxshift(yj-yi,boxysize)
+ zja=boxshift(zj-zi,boxzsize)
+ dist_init=xja**2+yja**2+zja**2
+ if (sqrt(dist_init).le.(10.0+r_buff_list)) then
+! Here the list is created
+! if (i.eq.2) then
+! print *,i,j,dist_init,ilist_catpnorm
+! endif
+ ilist_catpnorm=ilist_catpnorm+1
+
+! this can be substituted by cantor and anti-cantor
+ contlistcatpnormi(ilist_catpnorm)=i
+ contlistcatpnormj(ilist_catpnorm)=j
+ endif
+! enddo
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "before MPIREDUCE",ilist_catsctran,ilist_catptran,&
+ ilist_catscnorm,ilist_catpnorm,ilist_catscang
+
+ do i=1,ilist_catpnorm
+ write (iout,*) i,contlistcatpnormi(i)
+ enddo
+
+
+#endif
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(ilist_catpnorm,g_ilist_catcatnorm,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+ call MPI_Gather(ilist_catpnorm,1,MPI_INTEGER,&
+ i_ilist_catpnorm,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_ilist_catpnorm(i-1)+displ(i-1)
+ enddo
+! write(iout,*) "before gather",displ(0),displ(1)
+ call MPI_Gatherv(contlistcatpnormi,ilist_catpnorm,MPI_INTEGER,&
+ newcontlistcatcatnormi,i_ilist_catpnorm,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(contlistcatpnormj,ilist_catpnorm,MPI_INTEGER,&
+ newcontlistcatcatnormj,i_ilist_catpnorm,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Bcast(g_ilist_catcatnorm,1,MPI_INT,king,FG_COMM,IERR)
+! write(iout,*) "before bcast",g_ilist_sc
+! call MPI_Bcast(g_ilist_sc,1,MPI_INT,king,FG_COMM)
+ call MPI_Bcast(newcontlistcatcatnormi,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR)
+ call MPI_Bcast(newcontlistcatcatnormj,g_ilist_catcatnorm,MPI_INT,king,FG_COMM,IERR)
+
+
+ else
+ g_ilist_catcatnorm=ilist_catpnorm
+ do i=1,ilist_catpnorm
+ newcontlistcatcatnormi(i)=contlistcatpnormi(i)
+ newcontlistcatcatnormj(i)=contlistcatpnormj(i)
+ enddo
+ endif
+ call int_bounds(g_ilist_catcatnorm,g_listcatcatnorm_start,g_listcatcatnorm_end)
+
+#ifdef DEBUG
+ write (iout,*) "after MPIREDUCE",g_ilist_catcatnorm
+
+ do i=1,g_ilist_catcatnorm
+ write (iout,*) i,newcontlistcatcatnormi(i),newcontlistcatcatnormj(i)
+ enddo
+#endif
+ write(iout,*),"END make_catcat"
+ return
+ end subroutine make_cat_cat_list
!-----------------------------------------------------------------------------
return
end function sq
+#ifdef LBFGS
+ double precision function funcgrad(x,g)
+ use MD_data, only: totT,usampl
+ implicit none
+ double precision energia(0:n_ene)
+ double precision x(nvar),g(nvar)
+ integer i
+ call var_to_geom(nvar,x)
+ call zerograd
+ call chainbuild
+ call etotal(energia(0))
+ call sum_gradient
+ funcgrad=energia(0)
+ call cart2intgrad(nvar,g)
+ if (usampl) then
+ do i=1,nres-3
+ gloc(i,icg)=gloc(i,icg)+dugamma(i)
+ enddo
+ do i=1,nres-2
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+ enddo
+ endif
+ do i=1,nvar
+ g(i)=g(i)+gloc(i,icg)
+ enddo
+ return
+ end function funcgrad
+ subroutine cart2intgrad(n,g)
+ integer n
+ double precision g(n)
+ double precision drt(3,3,nres),rdt(3,3,nres),dp(3,3),&
+ temp(3,3),prordt(3,3,nres),prodrt(3,3,nres)
+ double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
+ double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,&
+ cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
+ double precision fromto(3,3),aux(6)
+ integer i,ii,j,jjj,k,l,m,indi,ind,ind1
+ logical sideonly
+ sideonly=.false.
+ g=0.0d0
+ if (sideonly) goto 10
+ do i=1,nres-2
+ rdt(1,1,i)=-rt(1,2,i)
+ rdt(1,2,i)= rt(1,1,i)
+ rdt(1,3,i)= 0.0d0
+ rdt(2,1,i)=-rt(2,2,i)
+ rdt(2,2,i)= rt(2,1,i)
+ rdt(2,3,i)= 0.0d0
+ rdt(3,1,i)=-rt(3,2,i)
+ rdt(3,2,i)= rt(3,1,i)
+ rdt(3,3,i)= 0.0d0
+ enddo
+ do i=2,nres-2
+ drt(1,1,i)= 0.0d0
+ drt(1,2,i)= 0.0d0
+ drt(1,3,i)= 0.0d0
+ drt(2,1,i)= rt(3,1,i)
+ drt(2,2,i)= rt(3,2,i)
+ drt(2,3,i)= rt(3,3,i)
+ drt(3,1,i)=-rt(2,1,i)
+ drt(3,2,i)=-rt(2,2,i)
+ drt(3,3,i)=-rt(2,3,i)
+ enddo
+ ind1=0
+ do i=1,nres-2
+ ind1=ind1+1
+ if (n.gt.nphi) then
+
+ do j=1,3
+ do k=1,2
+ dpjk=0.0D0
+ do l=1,3
+ dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
+ enddo
+ dp(j,k)=dpjk
+ prordt(j,k,i)=dp(j,k)
+ enddo
+ dp(j,3)=0.0D0
+ g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
+ enddo
+ xx1(1)=-0.5D0*xloc(2,i+1)
+ xx1(2)= 0.5D0*xloc(1,i+1)
+ do j=1,3
+ xj=0.0D0
+ do k=1,2
+ xj=xj+r(j,k,i)*xx1(k)
+ enddo
+ xx(j)=xj
+ enddo
+ do j=1,3
+ rj=0.0D0
+ do k=1,3
+ rj=rj+prod(j,k,i)*xx(k)
+ enddo
+ g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg)
+ enddo
+ if (i.lt.nres-2) then
+ do j=1,3
+ dxoiij=0.0D0
+ do k=1,3
+ dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+ enddo
+ g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg)
+ enddo
+ endif
+
+ endif
+
+
+ if (i.gt.1) then
+ do j=1,3
+ do k=1,3
+ dpjk=0.0
+ do l=2,3
+ dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
+ enddo
+ dp(j,k)=dpjk
+ prodrt(j,k,i)=dp(j,k)
+ enddo
+ g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
+ enddo
+ endif
+ xx(1)= 0.0D0
+ xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
+ xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
+ if (i.gt.1) then
+ do j=1,3
+ rj=0.0D0
+ do k=2,3
+ rj=rj+prod(j,k,i)*xx(k)
+ enddo
+ g(i-1)=g(i-1)-rj*gradx(j,i+1,icg)
+ enddo
+ endif
+ if (i.gt.1) then
+ do j=1,3
+ dxoiij=0.0D0
+ do k=1,3
+ dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+ enddo
+ g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg)
+ enddo
+ endif
+ do j=i+1,nres-2
+ ind1=ind1+1
+ call build_fromto(i+1,j+1,fromto)
+ do k=1,3
+ do l=1,3
+ tempkl=0.0D0
+ do m=1,2
+ tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
+ enddo
+ temp(k,l)=tempkl
+ enddo
+ enddo
+ if (n.gt.nphi) then
+ do k=1,3
+ g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
+ enddo
+ do k=1,3
+ dxoijk=0.0D0
+ do l=1,3
+ dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+ enddo
+ g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg)
+ enddo
+ endif
+ do k=1,3
+ do l=1,3
+ tempkl=0.0D0
+ do m=1,3
+ tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
+ enddo
+ temp(k,l)=tempkl
+ enddo
+ enddo
+ if (i.gt.1) then
+ do k=1,3
+ g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
+ enddo
+ do k=1,3
+ dxoijk=0.0D0
+ do l=1,3
+ dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+ enddo
+ g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg)
+ enddo
+ endif
+ enddo
+ enddo
+
+ if (nvar.le.nphi+ntheta) return
+
+ 10 continue
+ do i=2,nres-1
+ if (iabs(itype(i,1)).eq.10 .or. itype(i,1).eq.ntyp1& !) cycle
+ .or. mask_side(i).eq.0 ) cycle
+ ii=ialph(i,1)
+ dsci=vbld(i+nres)
+#ifdef OSF
+ alphi=alph(i)
+ omegi=omeg(i)
+ if(alphi.ne.alphi) alphi=100.0
+ if(omegi.ne.omegi) omegi=-100.0
+#else
+ alphi=alph(i)
+ omegi=omeg(i)
+#endif
+ cosalphi=dcos(alphi)
+ sinalphi=dsin(alphi)
+ cosomegi=dcos(omegi)
+ sinomegi=dsin(omegi)
+ temp(1,1)=-dsci*sinalphi
+ temp(2,1)= dsci*cosalphi*cosomegi
+ temp(3,1)=-dsci*cosalphi*sinomegi
+ temp(1,2)=0.0D0
+ temp(2,2)=-dsci*sinalphi*sinomegi
+ temp(3,2)=-dsci*sinalphi*cosomegi
+ theta2=pi-0.5D0*theta(i+1)
+ cost2=dcos(theta2)
+ sint2=dsin(theta2)
+ jjj=0
+ do j=1,2
+ xp=temp(1,j)
+ yp=temp(2,j)
+ xxp= xp*cost2+yp*sint2
+ yyp=-xp*sint2+yp*cost2
+ zzp=temp(3,j)
+ xx(1)=xxp
+ xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
+ xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
+ do k=1,3
+ dj=0.0D0
+ do l=1,3
+ dj=dj+prod(k,l,i-1)*xx(l)
+ enddo
+ aux(jjj+k)=dj
+ enddo
+ jjj=jjj+3
+ enddo
+ do k=1,3
+ g(ii)=g(ii)+aux(k)*gradx(k,i,icg)
+ g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg)
+ enddo
+ enddo
+ return
+ end subroutine cart2intgrad
+
+
+#endif
!--------------------------------------------------------------------------
end module energy
enddo
enddo
- call contact(.true.,ncont_ref,icont_ref,co)
+ call contact(.false.,ncont_ref,icont_ref,co)
! do k=1,nres
! write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
cref(j,i,kkk)=c(j,i)
enddo
enddo
- call contact(.true.,ncont_ref,icont_ref,co)
+ call contact(.false.,ncont_ref,icont_ref,co)
endif
! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
! call flush(iout)
if (pdbref) then
if(me.eq.king.or..not.out1file) &
write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
+ if (.false.) then
do i=1,ncont_ref
do j=1,2
icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
enddo
endif
+ endif
if (constr_homology.gt.0) then
! write (iout,*) "Calling read_constr_homology"
! call flush(iout)
write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
enddo
return
- end
+ end subroutine
!c---------------------------------------------------------------------
integer function tperm(i,iperm,tabpermchain)
! implicit none
!-----------------------------------------------------------------------------
subroutine read_x(kanal,*)
- use geometry_data
- use energy_data
- use geometry, only:int_from_cart1
+ use geometry_data
+ use energy_data
+ use geometry, only:int_from_cart1
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.GEO'
use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
use energy_data
-! use energy, only: to_box
! use control
use compare_data
use MD_data
ires=0
write(iout,*) "TUTUT"
do i=nnt,nct
- write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
+! write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
iti=itype(i,molnum(i))
print *,i,molnum(i)
if (molnum(i+1).eq.0) then
! if (zi.lt.0.0d0) zi=zi+boxzsize
write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
ires,xi,yi,zi,vtot(i)
- else
+ elseif (molnum(i).eq.5) then
+ xi=c(1,i)
+ yi=c(2,i)
+ zi=c(3,i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0.0d0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0.0d0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0.0d0) zi=zi+boxzsize
+ write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
+ ires,xi,yi,zi,vtot(i)
+ else
+
write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
ires,(c(j,i),j=1,3),vtot(i)
endif
enddo
write (iunit,'(a)') 'TER'
do i=nnt,nct-1
+ if (molnum(i).eq.5) cycle
if (itype(i,1).eq.ntyp1) cycle
if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
write (iunit,30) ica(i),ica(i+1)
!-----------------------------------------------------------------------------
integer function rescode(iseq,nam,itype,molecule)
-! use io_base, only: ucase
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.NAMES'
real(kind=8) :: energia(0:n_ene)
character(len=5) :: angid(4)=reshape((/'PHI ','THETA','ALPHA',&
'OMEGA'/),shape(angid))
- real(kind=8) :: ang_list(10)
+ real(kind=8) :: ang_list(10),ff
!el real(kind=8),dimension(6*nres) :: g,x !(maxvar) (maxvar=6*maxres)
real(kind=8),dimension(:),allocatable :: g,x !(maxvar) (maxvar=6*maxres)
integer :: nn(10)
!el local variables
- integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode,nfun
+ integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: etot,gnorm!,fdum
integer,dimension(1) :: uiparm
real(kind=8),dimension(1) :: urparm
etot = energia(0)
nf=1
nfl=3
- call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+#ifdef LBFGS
+ ff=funcgrad(x,g)
+#else
+
+!d write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
+ call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+!d write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar+20) !sp
+#endif
gnorm=0.0d0
do k=1,nvar
gnorm=gnorm+g(k)**2
! use csa_data
! use energy
implicit none
+!#ifdef LBFGS
+! double precision funcgrad_dc
+!#endif
!-----------------------------------------------------------------------------
!
!
! *** assess candidate step (***sol version 2.3) ***
!
integer :: liv, l
- integer(kind=8) ::lv
+ integer ::lv
integer :: iv(liv)
real(kind=8) :: v(lv)
!
! *** alg = 2 means general unconstrained optimization constants.
!
integer :: liv, l
- integer(kind=8)::lv
+ integer::lv
integer :: iv(liv)
integer :: alg
real(kind=8) :: v(lv)
! *** parameter declarations ***
!
integer :: liv, p
- integer(kind=8)::lv
+ integer::lv
integer :: iv(liv)
real(kind=8) :: d(p), g(p), v(lv), x(p)
!
! *** alg = 1 for regression, alg = 2 for general unconstrained opt.
!
integer :: alg, liv, n
- integer(kind=8) :: lv
+ integer :: lv
integer :: iv(liv)
real(kind=8) :: d(n), v(lv)
!
! *** alg = 2 means general unconstrained optimization constants.
!
integer :: alg, l
- integer(kind=8)::lv
+ integer::lv
real(kind=8) :: v(lv)
!/+
!el real(kind=8) :: dmax1
! *** (analytic) gradient and hessian provided by the caller. ***
!
integer :: liv, n
- integer(kind=8):: lv
+ integer:: lv
integer:: iv(liv)
integer :: uiparm(1)
real(kind=8) :: d(n), x(n), v(lv), urparm(1)
! *** parameter declarations ***
!
integer :: lh, liv, n
- integer(kind=8)::lv
+ integer::lv
integer :: iv(liv)
real(kind=8) :: d(n), fx, g(n), h(lh), v(lv), x(n)
!
! *** parameter declarations ***
!
integer :: liv, n
- integer(kind=8) ::lv
+ integer ::lv
integer :: iv(liv)
real(kind=8) :: d(n), hdiag(n), v(lv)
!
! minimize_p.F
!-----------------------------------------------------------------------------
subroutine minimize(etot,x,iretcode,nfun)
-
+#ifdef LBFGS
+ use minima
+ use inform
+ use output
+ use iounit
+ use scales
+ use energy, only: funcgrad,etotal,enerprint
+#else
use energy, only: func,gradient,fdum!,etotal,enerprint
+#endif
use comm_srutu
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
!el real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
- real(kind=8) :: energia(0:n_ene)
+ real(kind=8) :: energia(0:n_ene),grdmin
! external func,gradient,fdum
! external func_restr,grad_restr
logical :: not_done,change,reduce
!el common /przechowalnia/ v
!el local variables
integer :: iretcode,nfun,nvar_restr,idum(1),j
- integer (kind=8):: lv
+ integer :: lv
real(kind=8) :: etot,rdum(1) !,fdum
-
+#ifdef LBFGS
+ external optsave
+#endif
lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
-
+#ifndef LBFGS
if (.not.allocated(v)) allocate(v(1:lv))
+#endif
+#ifdef LBFGS
+ maxiter=maxmin
+ coordtype='RIGIDBODY'
+ grdmin=tolf
+ jout=iout
+ jprint=print_min_stat
+ 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
+ set_scale = .true.
+!c nvar = 0
+ do i = 1, nvar
+!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"
+ write (iout,*) 'Calling LBFGS, minimization in angles'
+ call var_to_geom(nvar,x)
+ call chainbuild
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ call lbfgs (nvar,x,etot,grdmin,funcgrad,optsave)
+ deallocate(scale)
+ write (iout,*) "Minimized energy",etot
+#else
icall = 1
! write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
! ENDDO ! NOT_DONE
-
+#endif
return
end subroutine minimize
!-----------------------------------------------------------------------------
! gradient_p.F
!-----------------------------------------------------------------------------
+#ifndef LBFGS
subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
use energy, only: cartder,zerograd,etotal,sum_gradient
!d enddo
return
end subroutine grad_restr
+#endif
+#ifndef LBFGS
!-----------------------------------------------------------------------------
subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
! endif
return
end subroutine func_restr
+#endif
!-----------------------------------------------------------------------------
! subroutine func(n,x,nf,f,uiparm,urparm,ufparm) in module energy
!-----------------------------------------------------------------------------
!el subroutine xx2x(x,xx) in module math
!-----------------------------------------------------------------------------
subroutine minim_dc(etot,iretcode,nfun)
+#ifdef LBFGS
+ use minima
+ use inform
+ use output
+ use iounit
+ use scales
+#endif
use MPI_data
- use energy, only: fdum,check_ecartint
+ use energy, only: fdum,check_ecartint,etotal,enerprint
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MPI
! include 'COMMON.MINIM'
! include 'COMMON.CHAIN'
integer :: iretcode,nfun,k,i,j,idum(1)
- integer(kind=8) :: lv
+ integer :: lv,nvarx,nf
integer,dimension(liv) :: iv
real(kind=8) :: minval !,v(1:77+(6*nres)*(6*nres+17)/2) !(1:lv)
- real(kind=8),dimension(6*nres) :: x,d,xx !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(:),allocatable :: x,d,xx !(maxvar) (maxvar=6*maxres)
!el common /przechowalnia/ v
- real(kind=8) :: energia(0:n_ene)
+ real(kind=8) :: energia(0:n_ene),grdmin
! external func_dc,grad_dc ,fdum
logical :: not_done,change,reduce
real(kind=8) :: g(6*nres),f1,etot,rdum(1) !,fdum
+#ifdef LBFGS
+ external optsave
+ maxiter=maxmin
+ coordtype='CARTESIAN'
+ grdmin=tolf
+ jout=iout
+ jprint=print_min_stat
+ iwrite=0
+ write(iout,*) "in minim_dc LBFGS"
+#else
lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
print *,"lv value",lv
do i=1,6*nres
d(i)=1.0D-1
enddo
-
+#endif
+ if (.not.allocated(x)) then
+ allocate(x(6*nres))
+ allocate(xx(6*nres))
+ allocate(d(6*nres))
+ endif
k=0
do i=1,nres-1
do j=1,3
enddo
endif
enddo
+ nvarx=k
+ call chainbuild_cart
+ call etotal(energia(0))
+ call enerprint(energia(0))
+#ifdef LBFGS
+!c
+!c From tinker
+!c
+!c perform dynamic allocation of some global arrays
+!c
+ if (.not. allocated(scale)) allocate (scale(nvarx))
+!c
+!c set scaling parameter for function and derivative values;
+!c use square root of median eigenvalue of typical Hessian
+!c
+ set_scale = .true.
+!c nvar = 0
+ do i = 1, nvarx
+!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
+ write (iout,*) "minim_dc Calling lbfgs"
+ call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
+ deallocate(scale)
+#else
! call check_ecartint
call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
! call check_ecartint
+#endif
k=0
do i=1,nres-1
do j=1,3
! (c(j,i+nres),j=1,3)
! enddo
!el----------------------------
+#ifndef LBFGS
etot=v(10)
iretcode=iv(1)
nfun=iv(6)
+#endif
return
end subroutine minim_dc
+
+
+#ifdef LBFGS
+ real(kind=8) function funcgrad_dc(x,g)
+ use energy, only: etotal, zerograd,cartgrad
+! implicit real*8 (a-h,o-z)
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ integer k,nf,i,j
+ real(kind=8),dimension(nres*6):: x,g
+ double precision energia(0:n_ene)
+! common /gacia/ nf
+!c
+ nf=nf+1
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ dc(j,i)=x(k)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ dc(j,i+nres)=x(k)
+ enddo
+ endif
+ enddo
+ call chainbuild_cart
+ call zerograd
+ call etotal(energia(0))
+!c write (iout,*) "energia",energia(0)
+ funcgrad_dc=energia(0)
+ call cartgrad
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ g(k)=gcart(j,i)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ g(k)=gxcart(j,i)
+ enddo
+ endif
+ enddo
+
+ return
+ end function
+#else
!-----------------------------------------------------------------------------
subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
return
end subroutine grad_dc
+#endif
!-----------------------------------------------------------------------------
! minim_mcmf.F
!-----------------------------------------------------------------------------
use MPI_data
use csa_data
+#ifndef LBFGS
use energy, only: func,gradient,fdum
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
integer,dimension(6) :: indx
integer,dimension(liv) :: iv
integer :: idum(1),nf !
- integer(kind=8) :: lv
+ integer :: lv
real(kind=8) :: rdum(1)
real(kind=8) :: przes(3),obrot(3,3),eee
logical :: non_conv
call mpi_send(ene0,1,mpi_double_precision,king,idreal,&
CG_COMM,ierr)
go to 10
+#endif
return
end subroutine minim_mcmf
#endif
end subroutine sc_minimize
!-----------------------------------------------------------------------------
subroutine minimize_sc1(etot,x,iretcode,nfun)
-
+#ifndef LBFGS
use energy, only: func,gradient,fdum,etotal,enerprint
+#endif
use comm_srutu
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! external func_restr1,grad_restr1
logical :: not_done,change,reduce
!el common /przechowalnia/ v
- integer(kind=8):: lv
+ integer:: lv
integer :: nvar_restr,i,j
integer :: idum(1)
real(kind=8) :: rdum(1),etot !,fdum
-
+#ifndef LBFGS
lv=(77+(6*nres)*(6*nres+17)/2) !(77+maxvar*(maxvar+17)/2)) (maxvar=6*maxres)
if (.not. allocated(v)) allocate(v(1:lv))
etot=v(10)
iretcode=iv(1)
nfun=iv(6)
-
+#endif
return
end subroutine minimize_sc1
!-----------------------------------------------------------------------------
return
end subroutine func_restr1
!-----------------------------------------------------------------------------
+
subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
use energy, only: cartder,zerograd,esc,sc_grad
!
! use control
integer :: n, liv
- integer(kind=8)::lv
+ integer::lv
integer:: iv(liv)
integer :: uiparm(1)
real(kind=8) :: d(n), x(n), v(lv), urparm(1)
!
integer :: iv1, nf
real(kind=8) :: f
- integer(kind=8)::g1
+ integer::g1
!
! *** subscripts for iv ***
!
! *** parameter declarations ***
!
integer :: liv, n
- integer(kind=8) :: lv
+ integer :: lv
integer :: iv(liv)
real(kind=8) :: d(n), fx, g(n), v(lv), x(n)
!
! *** parameter declarations ***
!
integer :: n
- integer(kind=8) :: lv
+ integer :: lv
real(kind=8) :: dig(n), nwtstp(n), step(n), v(lv)
!
! *** purpose ***
! *** last card of wzbfgs follows ***
end subroutine wzbfgs
!-----------------------------------------------------------------------------
+
+!
+!
+! ###################################################
+! ## COPYRIGHT (C) 1999 by Jay William Ponder ##
+! ## All Rights Reserved ##
+! ###################################################
+!
+! ##############################################################
+! ## ##
+! ## subroutine lbfgs -- limited memory BFGS optimization ##
+! ## ##
+! ##############################################################
+!
+!
+! "lbfgs" is a limited memory BFGS quasi-newton nonlinear
+! optimization routine
+!
+! literature references:
+!
+! J. Nocedal, "Updating Quasi-Newton Matrices with Limited
+! Storage", Mathematics of Computation, 35, 773-782 (1980)
+!
+! D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
+! for Large Scale Optimization", Mathematical Programming,
+! 45, 503-528 (1989)
+!
+! J. Nocedal and S. J. Wright, "Numerical Optimization",
+! Springer-Verlag, New York, 1999, Section 9.1
+!
+! variables and parameters:
+!
+! nvar number of parameters in the objective function
+! x0 contains starting point upon input, upon return
+! contains the best point found
+! minimum during optimization contains best current function
+! value; returns final best function value
+! grdmin normal exit if rms gradient gets below this value
+! ncalls total number of function/gradient evaluations
+!
+! required external routines:
+!
+! fgvalue function to evaluate function and gradient values
+! optsave subroutine to write out info about current status
+!
+!
+ subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
+ use inform
+ use iounit
+ use keys
+ use linmin
+ use math2
+ use minima
+ use output
+ use scales
+ use MD_data, only: itime_mat
+#ifdef LBFGS
+ use energy_data,only:statusbf,niter,ncalls
+#endif
+ implicit none
+ integer i,j,k,m
+ integer nvar,next
+ integer msav,muse
+#ifndef LBFGS
+ integer niter,ncalls
+ character*9 statusbf
+#endif
+ integer nerr,maxerr
+ real*8 f,f_old,fgvalue
+ real*8 f_move,x_move
+ real*8 g_norm,g_rms
+ real*8 minimum,grdmin
+ real*8 angle,rms,beta
+ real*8 ys,yy,gamma
+ real*8 x0(*)
+ real*8, allocatable :: rho(:)
+ real*8, allocatable :: alpha(:)
+ real*8, allocatable :: x_old(:)
+ real*8, allocatable :: g(:)
+ real*8, allocatable :: g_old(:)
+ real*8, allocatable :: p(:)
+ real*8, allocatable :: q(:)
+ real*8, allocatable :: r(:)
+ real*8, allocatable :: h0(:)
+ real*8, allocatable :: s(:,:)
+ real*8, allocatable :: y(:,:)
+ logical done
+ character*9 blank
+ character*20 keyword
+ character*240 record
+ character*240 string
+ external fgvalue,optsave
+! common /lbfgstat/ status,niter,ncalls
+!c
+!c
+!c initialize some values to be used below
+!c
+ ncalls = 0
+ rms = sqrt(dble(nvar))
+ if (coordtype .eq. 'CARTESIAN') then
+ rms = rms / sqrt(3.0d0)
+ else if (coordtype .eq. 'RIGIDBODY') then
+ rms = rms / sqrt(6.0d0)
+ end if
+ blank = ' '
+ done = .false.
+ nerr = 0
+ maxerr = 2
+!c
+!c perform dynamic allocation of some global arrays
+!c
+ if (.not. allocated(scale)) allocate (scale(nvar))
+!c
+!c set default values for variable scale factors
+!c
+ if (.not. set_scale) then
+ do i = 1, nvar
+ if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0
+ end do
+ end if
+!c
+!c set default parameters for the optimization
+!c
+ msav = min(nvar,20)
+ if (fctmin .eq. 0.0d0) fctmin = -100000000.0d0
+ if (maxiter .eq. 0) maxiter = 1000000
+ if (nextiter .eq. 0) nextiter = 1
+ if (jprint .lt. 0) jprint = 1
+ if (iwrite .lt. 0) iwrite = 1
+ write(iout,*) "maxiter",maxiter
+!c
+!c set default parameters for the line search
+!c
+ if (stpmax .eq. 0.0d0) stpmax = 5.0d0
+ stpmin = 1.0d-16
+ cappa = 0.9d0
+ slpmax = 10000.0d0
+ angmax = 180.0d0
+ intmax = 5
+!c
+!c search the keywords for optimization parameters
+!c
+#ifdef LBFGSREAD
+ do i = 1, nkey
+ next = 1
+ record = keyline(i)
+ call gettext (record,keyword,next)
+ call upcase (keyword)
+ string = record(next:240)
+ if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
+ read (string,*,err=10,end=10) msav
+ msav = max(0,min(msav,nvar))
+ else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
+ msav = 0
+ else if (keyword(1:7) .eq. 'FCTMIN ') then
+ read (string,*,err=10,end=10) fctmin
+ else if (keyword(1:8) .eq. 'MAXITER ') then
+ read (string,*,err=10,end=10) maxiter
+ else if (keyword(1:8) .eq. 'STEPMAX ') then
+ read (string,*,err=10,end=10) stpmax
+ else if (keyword(1:8) .eq. 'STEPMIN ') then
+ read (string,*,err=10,end=10) stpmin
+ else if (keyword(1:6) .eq. 'CAPPA ') then
+ read (string,*,err=10,end=10) cappa
+ else if (keyword(1:9) .eq. 'SLOPEMAX ') then
+ read (string,*,err=10,end=10) slpmax
+ else if (keyword(1:7) .eq. 'ANGMAX ') then
+ read (string,*,err=10,end=10) angmax
+ else if (keyword(1:7) .eq. 'INTMAX ') then
+ read (string,*,err=10,end=10) intmax
+ end if
+ 10 continue
+ end do
+#endif
+!c
+!c print header information about the optimization method
+!c
+ if (jprint .gt. 0) then
+ if (msav .eq. 0) then
+ write (jout,20)
+ 20 format (/,' Steepest Descent Gradient Optimization :')
+ write (jout,30)
+ 30 format (/,' SD Iter F Value G RMS F Move',&
+ ' X Move Angle FG Call Comment',/)
+ else
+ write (jout,40)
+ 40 format (/,' Limited Memory BFGS Quasi-Newton',&
+ ' Optimization :')
+ write (jout,50)
+ 50 format (/,' QN Iter F Value G RMS F Move',&
+ ' X Move Angle FG Call Comment',/)
+ end if
+ flush (jout)
+ end if
+!c
+!c perform dynamic allocation of some local arrays
+!c
+ allocate (x_old(nvar))
+ allocate (g(nvar))
+ allocate (g_old(nvar))
+ allocate (p(nvar))
+ allocate (q(nvar))
+ allocate (r(nvar))
+ allocate (h0(nvar))
+ if (msav .ne. 0) then
+ allocate (rho(msav))
+ allocate (alpha(msav))
+ allocate (s(nvar,msav))
+ allocate (y(nvar,msav))
+ end if
+!c
+!c evaluate the function and get the initial gradient
+!c
+ niter = nextiter - 1
+ maxiter = niter + maxiter
+ ncalls = ncalls + 1
+! itime_mat=niter
+ f = fgvalue (x0,g)
+ f_old = f
+ m = 0
+ gamma = 1.0d0
+ g_norm = 0.0d0
+ g_rms = 0.0d0
+ do i = 1, nvar
+ g_norm = g_norm + g(i)*g(i)
+ g_rms = g_rms + (g(i)*scale(i))**2
+ end do
+ g_norm = sqrt(g_norm)
+ g_rms = sqrt(g_rms) / rms
+ f_move = 0.5d0 * stpmax * g_norm
+!c
+!c print initial information prior to first iteration
+!c
+ write(jout,*) "before first"
+ if (jprint .gt. 0) then
+ if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
+ write (jout,60) niter,f,g_rms,ncalls
+ 60 format (i6,f14.4,f11.4,2x,i7)
+ else
+ write (jout,70) niter,f,g_rms,ncalls
+ 70 format (i6,d14.4,d11.4,2x,i7)
+ end if
+ flush (jout)
+ end if
+!c
+!c write initial intermediate prior to first iteration
+!c
+ if (iwrite .gt. 0) call optsave (niter,f,x0)
+!c
+!c tests of the various termination criteria
+!c
+ if (niter .ge. maxiter) then
+ statusbf = 'IterLimit'
+ done = .true.
+ end if
+ if (f .le. fctmin) then
+ statusbf = 'SmallFct '
+ done = .true.
+ end if
+ if (g_rms .le. grdmin) then
+ statusbf = 'SmallGrad'
+ done = .true.
+ end if
+!c
+!c start of a new limited memory BFGS iteration
+!c
+ do while (.not. done)
+ niter = niter + 1
+!c write (jout,*) "LBFGS niter",niter
+ muse = min(niter-1,msav)
+ m = m + 1
+ if (m .gt. msav) m = 1
+!c
+!c estimate Hessian diagonal and compute the Hg product
+!c
+ do i = 1, nvar
+ h0(i) = gamma
+ q(i) = g(i)
+ end do
+ k = m
+ do j = 1, muse
+ k = k - 1
+ if (k .eq. 0) k = msav
+ alpha(k) = 0.0d0
+ do i = 1, nvar
+ alpha(k) = alpha(k) + s(i,k)*q(i)
+ end do
+ alpha(k) = alpha(k) * rho(k)
+ do i = 1, nvar
+ q(i) = q(i) - alpha(k)*y(i,k)
+ end do
+ end do
+ do i = 1, nvar
+ r(i) = h0(i) * q(i)
+ end do
+ do j = 1, muse
+ beta = 0.0d0
+ do i = 1, nvar
+ beta = beta + y(i,k)*r(i)
+ end do
+ beta = beta * rho(k)
+ do i = 1, nvar
+ r(i) = r(i) + s(i,k)*(alpha(k)-beta)
+ end do
+ k = k + 1
+ if (k .gt. msav) k = 1
+ end do
+!c
+!c set search direction and store current point and gradient
+!c
+ do i = 1, nvar
+ p(i) = -r(i)
+ x_old(i) = x0(i)
+ g_old(i) = g(i)
+ end do
+!c
+!c perform line search along the new conjugate direction
+!c
+ statusbf = blank
+!c write (jout,*) "Before search"
+ call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,statusbf)
+!c write (jout,*) "After search"
+!c
+!c update variables based on results of this iteration
+!c
+ if (msav .ne. 0) then
+ ys = 0.0d0
+ yy = 0.0d0
+ do i = 1, nvar
+ s(i,m) = x0(i) - x_old(i)
+ y(i,m) = g(i) - g_old(i)
+ ys = ys + y(i,m)*s(i,m)
+ yy = yy + y(i,m)*y(i,m)
+ end do
+ gamma = abs(ys/yy)
+ rho(m) = 1.0d0 / ys
+ end if
+!c
+!c get the sizes of the moves made during this iteration
+!c
+ f_move = f_old - f
+! if (f_move.eq.0) f_move=1.0d0
+ f_old = f
+ x_move = 0.0d0
+ do i = 1, nvar
+ x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
+ end do
+ x_move = sqrt(x_move) / rms
+ if (coordtype .eq. 'INTERNAL') then
+ x_move = radian * x_move
+ end if
+!c
+!c compute the rms gradient per optimization parameter
+!c
+ g_rms = 0.0d0
+ do i = 1, nvar
+ g_rms = g_rms + (g(i)*scale(i))**2
+ end do
+ g_rms = sqrt(g_rms) / rms
+!c
+!c test for error due to line search problems
+!c
+ if (statusbf.eq.'BadIntpln' .or. statusbf.eq.'IntplnErr') then
+ nerr = nerr + 1
+ if (nerr .ge. maxerr) done = .true.
+ else
+ nerr = 0
+ end if
+!c
+!c test for too many total iterations
+!c
+ if (niter .ge. maxiter) then
+ statusbf = 'IterLimit'
+ done = .true.
+ end if
+!c
+!c test the normal termination criteria
+!c
+ if (f .le. fctmin) then
+ statusbf = 'SmallFct '
+ done = .true.
+ end if
+ if (g_rms .le. grdmin) then
+ statusbf = 'SmallGrad'
+ done = .true.
+ end if
+!c
+!c print intermediate results for the current iteration
+!c
+! write(iout,*) "niter1",niter
+ itime_mat=niter
+ if (jprint .gt. 0) then
+ if (done .or. mod(niter,jprint).eq.0) then
+ if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.&
+ g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.&
+ f_move.gt.-1.0d5) then
+ write (jout,80) niter,f,g_rms,f_move,x_move,&
+ angle,ncalls,statusbf
+ 80 format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
+ else
+ write (jout,90) niter,f,g_rms,f_move,x_move,&
+ angle,ncalls,statusbf
+ 90 format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
+ end if
+ end if
+ flush (jout)
+ end if
+!c
+!c write intermediate results for the current iteration
+!c
+ if (iwrite .gt. 0) then
+ if (done .or. mod(niter,iwrite).eq.0) then
+ call optsave (niter,f,x0)
+ end if
+ end if
+ end do
+!c
+!c perform deallocation of some local arrays
+!c
+ deallocate (x_old)
+ deallocate (g)
+ deallocate (g_old)
+ deallocate (p)
+ deallocate (q)
+ deallocate (r)
+ deallocate (h0)
+ if (msav .ne. 0) then
+ deallocate (rho)
+ deallocate (alpha)
+ deallocate (s)
+ deallocate (y)
+ end if
+!c
+!c set final value of the objective function
+!c
+ minimum = f
+ if (jprint .gt. 0) then
+ if (statusbf.eq.'SmallGrad' .or. statusbf.eq.'SmallFct ') then
+ write (jout,100) statusbf
+ 100 format (/,' LBFGS -- Normal Termination due to ',a9)
+ else
+ write (jout,110) statusbf
+ 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)
+ end if
+ flush (jout)
+ end if
+ return
+ end subroutine
+!ifdef LBFGS
+! double precision function funcgrad(x,g)
+! implicit none
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.MD'
+! include 'COMMON.QRESTR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! double precision energia(0:n_ene)
+! double precision x(nvar),g(nvar)
+! integer i
+! if (jjj.gt.0) then
+! write (iout,*) "in func x"
+! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+! endif
+! call var_to_geom(nvar,x)
+! call zerograd
+! call chainbuild_extconf
+! call etotal(energia(0))
+! call sum_gradient
+! funcgrad=energia(0)
+! call cart2intgrad(nvar,g)
+!
+! Add the components corresponding to local energy terms.
+!
+! Add the usampl contributions
+! if (usampl) then
+! do i=1,nres-3
+! gloc(i,icg)=gloc(i,icg)+dugamma(i)
+! enddo
+! do i=1,nres-2
+! gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+! enddo
+! endif
+! do i=1,nvar
+!d write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
+! g(i)=g(i)+gloc(i,icg)
+! enddo
+! return
+! end
+!endif
!-----------------------------------------------------------------------------
end module minimm
+
call exec_softreg
else if (modecalc.eq.12) then
call exec_MD
- call exec_checkgrad
+! call exec_checkgrad
else if (modecalc.eq.14) then
call exec_MREMD
else
real(kind=8) :: energy_long(0:n_ene),energy_short(0:n_ene)
real(kind=8) :: varia(6*nres) !(maxvar) (maxvar=6*maxres)
real(kind=8) :: time00, evals, etota, etot, time_ene, time1
- integer :: nharp,nft_sc,iretcode,nfun
+ integer :: nharp,nft_sc,iretcode
+#ifndef LBFGS
+ integer :: nfun
+#endif
integer,dimension(4,nres) :: iharp !(4,nres/3)(4,maxres/3)
logical :: fail
real(kind=8) :: rms,frac,frac_nn,co
! call check_eint
! call exec_checkgrad !el
endif
+! print *,'SUMSL return code is',iretcode,' eval ',nfun
+#ifdef LBFGS
+ print *,'LBFGS return code is',status,' eval ',nfun
+#else
print *,'SUMSL return code is',iretcode,' eval ',nfun
+#endif
+
#ifdef MPI
evals=nfun/(MPI_WTIME()-time1)
#endif
use geometry_data !include 'COMMON.GEO''COMMON.CHAIN'
! use REMD !include 'COMMON.REMD'
! use MD !include 'COMMON.MD'
+
use regularize_
use compare
implicit none
if (outpdb) call pdbout(etot,titel(:32),ipdb)
if (outmol2) call mol2out(etot,titel(:32))
if (refstr) call rms_nac_nnc(rms,frac,frac_nn,co,.true.)
+#ifdef LBFGS
+ write (iout,'(a,a9)') 'LBFGS return code:',status
+#else
write (iout,'(a,i3)') 'SUMSL return code:',iretcode
+#endif
return
end subroutine exec_regularize
!-----------------------------------------------------------------------------
real(kind=8) :: energy_(0:n_ene)
logical :: eof
real(kind=8) :: etot,ene0
- integer :: mm,imm,nft,n,iconf,nmin,i,iretcode,nfun,it,&
+ integer :: mm,imm,nft,n,iconf,nmin,i,iretcode,it,&
nf_mcmf,j
+#ifndef LBFGS
+ integer :: nfun
+#endif
real(kind=8) :: rms,frac,frac_nn,co,time,ene
eof=.false.
projects / unres4.git / commitdiff