d_t_new(j,0)=d_t_old(j,0)+adt2
d_t(j,0)=d_t_old(j,0)+adt
enddo
- do i=nnt,nct-1
+ do i=nnt,nct-1
C SPYTAC ADAMA
C do i=0,nres
+#ifdef DEBUG
+ write (iout,*) "i",i," d_a_old",(d_a_old(j,i),j=1,3)
+#endif
do j=1,3
adt=d_a_old(j,i)*d_time
adt2=0.5d0*adt
include 'COMMON.NAMES'
include 'COMMON.REMD'
include 'COMMON.TIME1'
+#ifdef LBFGS
+ character*9 status
+ integer niter
+ common /lbfgstat/ status,niter,nfun
+#endif
+ integer n_model_try,list_model_try(max_template),k
double precision tt0
real*8 energia_long(0:n_ene),
& energia_short(0:n_ene),vcm(3),incr(3)
endif
call flush(iout)
write (iout,*) "init_MD before initial structure REST ",rest
- if (.not.rest) then
+ if (.not.rest) then
+ 122 continue
if (iranconf.ne.0) then
c 8/22/17 AL Loop to produce a low-energy random conformation
do iranmin=1,10
44 continue
else if (preminim) then
if (start_from_model) then
- i_model=iran_num(1,constr_homology)
- write (iout,*) 'starting from model ',i_model
- do i=1,2*nres
- do j=1,3
- c(j,i)=chomo(j,i,i_model)
+ n_model_try=0
+ do while (fail .and. n_model_try.lt.constr_homology)
+ do
+ i_model=iran_num(1,constr_homology)
+ do k=1,n_model_try
+ if (i_model.eq.list_model_try(k)) exit
+ enddo
+ if (k.gt.n_model_try) exit
enddo
- enddo
- call int_from_cart(.true.,.false.)
- call sc_loc_geom(.false.)
- do i=1,nres-1
- do j=1,3
- dc(j,i)=c(j,i+1)-c(j,i)
- dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ n_model_try=n_model_try+1
+ list_model_try(n_model_try)=i_model
+ write (iout,*) 'starting from model ',i_model
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=chomo(j,i,i_model)
+ enddo
enddo
- enddo
- do i=2,nres-1
- do j=1,3
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
enddo
- enddo
- endif
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+ enddo
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies before removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
! Remove SC overlaps if requested
- if (overlapsc) then
- write (iout,*) 'Calling OVERLAP_SC'
- call overlap_sc(fail)
- endif
+ if (overlapsc) then
+ write (iout,*) 'Calling OVERLAP_SC'
+ call overlap_sc(fail)
+ if (fail) then
+ write (iout,*)
+ & "Failed to remove overlap from model",i_model
+ cycle
+ endif
+ endif
+#ifdef SEARCHSC
+ if (me.eq.king.or..not.out1file) then
+ write (iout,*) "Energies after removing overlaps"
+ call etotal(energia(0))
+ call enerprint(energia(0))
+ endif
! Search for better SC rotamers if requested
- if (searchsc) then
- call sc_move(2,nres-1,10,1d10,nft_sc,etot)
- print *,'SC_move',nft_sc,etot
- if (me.eq.king.or..not.out1file)
- & write(iout,*) 'SC_move',nft_sc,etot
+ if (searchsc) then
+ call sc_move(2,nres-1,10,1d10,nft_sc,etot)
+ print *,'SC_move',nft_sc,etot
+ if (me.eq.king.or..not.out1file)
+ & write(iout,*) 'SC_move',nft_sc,etot
+ endif
+ call etotal(energia(0))
+#endif
+ enddo
+ if (n_model_try.gt.constr_homology) then
+ write (iout,*)
+ & "All models have irreparable overlaps. Trying randoms starts."
+ iranconf=1
+ goto 122
+ endif
endif
- call etotal(energia(0))
C 8/22/17 AL Minimize initial structure
if (dccart) then
if (me.eq.king.or..not.out1file) write(iout,*)
- & 'Minimizing initial PDB structure: Calling MINIM_DC'
+ & 'Minimizing initial PDB structure: Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
else
call geom_to_var(nvar,varia)
if(me.eq.king.or..not.out1file) write (iout,*)
- & 'Minimizing initial PDB structure: Calling MINIMIZE.'
+ & 'Minimizing initial PDB structure: Calling MINIMIZE.'
call minimize(etot,varia,iretcode,nfun)
call var_to_geom(nvar,varia)
- endif
- if (me.eq.king.or..not.out1file)
+#ifdef LBFGS
+ if (me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+ if(me.eq.king.or..not.out1file)
+ & write(iout,*) 'LBFGS return code is ',status,' eval ',nfun
+#else
+ if (me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
- if(me.eq.king.or..not.out1file)
+ if(me.eq.king.or..not.out1file)
& write(iout,*) 'SUMSL return code is',iretcode,' eval ',nfun
+#endif
+ endif
endif
- endif
+ endif ! .not. rest
call chainbuild_cart
call kinetic(EK)
if (tbf) then
call verlet_bath
- endif
+ endif
kinetic_T=2.0d0/(dimen3*Rb)*EK
if(me.eq.king.or..not.out1file)then
call cartprint
& d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
d_t(:,chain_border(2,ichain-1))=0.0d0
enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+ & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
#else
ibeg=0
c do j=1,3
projects / unres.git / blobdiff