set (CMAKE_Fortran_FLAGS_RELEASE " ")
set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g -traceback")
# set(FFLAGS0 "-fpp -c -CB -g -ip " )
- set(FFLAGS0 "-CB -g -ip -fpp -mcmodel=large" )
+ set(FFLAGS0 "-O3 -ip -fpp -mcmodel=large" )
# set(FFLAGS0 "-O0 -CB -CA -g" )
set(FFLAGS1 "-fpp -c -O " )
set(FFLAGS2 "-fpp -c -g -CA -CB ")
integer :: i,j,k,iti,mnum,term
real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
mag1,mag2,v(3)
-#ifdef DEBUG
+!#ifdef DEBUG
write (iout,*) "Velocities, kietic"
do i=0,nres
write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),&
(d_t(j,i+nres),j=1,3)
enddo
-#endif
+!#endif
KEt_p=0.0d0
KEt_sc=0.0d0
! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres)
do i=nnt,term
mnum=molnum(i)
if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
-! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
+ write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3),mp(mnum)
if (mnum.gt.4) then
do j=1,3
v(j)=incr(j)+0.5d0*d_t(j,i)
incr(j)=incr(j)+d_t(j,i)
enddo
enddo
-! write(iout,*) 'KEt_p', KEt_p
+ write(iout,*) 'KEt_p', KEt_p
! The translational part for the side chain virtual bond
! Only now we can initialize incr with zeros. It must be equal
! to the velocities of the first Calpha.
v(j)=incr(j)+d_t(j,nres+i)
enddo
endif
-! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
+ write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3)
! write (iout,*) "i",i," msc",msc(iti,mnum)," v",(v(j),j=1,3)
KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
enddo
enddo
! goto 111
-! write(iout,*) 'KEt_sc', KEt_sc
+ write(iout,*) 'KEt_sc', KEt_sc
! The part due to stretching and rotation of the peptide groups
KEr_p=0.0D0
do i=nnt,nct-1
do j=1,3
incr(j)=d_t(j,i)
enddo
-! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3)
+ write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3),KEr_p,Ip(mnum)
KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) &
+incr(3)*incr(3))
enddo
! goto 111
-! write(iout,*) 'KEr_p', KEr_p
+ write(iout,*) 'KEr_p', KEr_p
! The rotational part of the side chain virtual bond
KEr_sc=0.0D0
do i=nnt,nct
do j=1,3
incr(j)=d_t(j,nres+i)
enddo
-! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
+ write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3)
KEr_sc=KEr_sc+Isc(iti,mnum)*(incr(1)*incr(1)+incr(2)*incr(2)+ &
incr(3)*incr(3))
endif
enddo
! The total kinetic energy
111 continue
-! write(iout,*) 'KEr_sc', KEr_sc
+ write(iout,*) 'KEr_sc', KEr_sc
KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
-! write (iout,*) "KE_total",KE_total
+ write (iout,*) "KE_total",KE_total
return
end subroutine kinetic
!-----------------------------------------------------------------------------
enddo ! ichain
!c The total kinetic energy
111 continue
-!c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,
-!c & ' KEr_sc', KEr_sc
+ write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p,&
+ ' KEr_sc', KEr_sc
KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc)
!c write (iout,*) "KE_total",KE_tota
#else
itime_mat=itime
if (ntwe.ne.0) then
if (mod(itime,ntwe).eq.0) then
-! call returnbox
+ call returnbox
call statout(itime)
! call returnbox
! call check_ecartint
! write (iout,*) "before enter if",scalel,icount_scale
if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
! write(iout,*) "I enter if"
+! energy_dec=.true.
+! call check_ecartint
! Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
scalel=.true.
ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax)) &
! call ginv_mult(fric_work, d_af_work)
! call ginv_mult(stochforcvec, d_as_work)
#ifdef FIVEDIAG
+#ifdef DEBUG
write(iout,*) "forces before fivediaginv"
do i=1,dimen*3
write(iout,*) "fricwork",i,fric_work(i)
enddo
+#endif
call fivediaginv_mult(dimen,fric_work, d_af_work)
call fivediaginv_mult(dimen,stochforcvec, d_as_work)
if (large) then
write(iout,*) "Initial state will be read from file ",&
rest2name(:ilen(rest2name))
call readrst
+ call chainbuild_cart
+ do i= 1,nres
+ do j=1,3
+ write(iout,*),"myc",c(j,i),dc(j,i)
+ enddo
+ enddo
endif
call rescale_weights(t_bath)
else
if(dccart)then
print *, 'Calling MINIM_DC'
call minim_dc(etot,iretcode,nfun)
+ print *,"fter MINIM_DC"
call int_from_cart1(.false.)
else
call geom_to_var(nvar,varia)
endif
enddo
endif
- endif
+ endif
+
call chainbuild_cart
call kinetic(EK)
write(iout,*) "just after kinetic"
ind=ind+3
do i=3,n
ind=ind+i-3
+#ifdef DEBUG
write (iout,*) "i",i," ind",ind," indu2",innt+i-2,&
" indu1",innt+i-1," indm",innt+i
+ write(iout,*) "value",innt-1+i-2,du2orig(innt-1+i-2),innt-1+i-1,du1orig(innt-1+i-1),innt-1+i,dmorig(innt-1+i)
+#endif
ghalf(ind+1)=du2orig(innt-1+i-2)
ghalf(ind+2)=du1orig(innt-1+i-1)
ghalf(ind+3)=dmorig(innt-1+i)
enddo
#endif
#ifdef DEBUG
+
write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
write (iout,*) "Five-diagonal inertia matrix, lower triangle"
-! call matoutr(n,ghalf)
+! call matout(n,ghalf)
+ write (iout,*), "MY IND",ind
+ do i=1,ind
+ write(iout,*) i,ghalf(i)
+ enddo
+#endif
+ mnum=molnum(innt+1)
+ call gldiag(1300*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
+ write(iout,*) "after internal",ierr,iwork
+#ifdef DEBUG
+ do i=1,n
+ do j=1,n
+ if (Gvec(i,j).ne.0) write(iout,*) i,j,Gvec(i,j),"kupa"
+ enddo
+ enddo
#endif
- call gldiag(nres*2,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
if (large) then
write (iout,'(//a,i3)')&
"Eigenvectors and eigenvalues of the G matrix chain",ichain
do i=1,n
do k=1,3
ii=ii+1
- mnum=molnum(innt_org)
- if (molnum(innt_org).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
+ mnum=molnum(innt_org+1)
+ if (molnum(innt_org+1).ge.4) geigen(i)=3.0/msc(itype(innt_org+i-1,mnum),mnum)
+!should not it be other way around
! if (molnum(innt).eq.5) write(iout,*) "typ",i,innt-1+i,itype(innt+i-1,5)
+
sigv=dsqrt((Rb*t_bath)/geigen(i))
lowb=-5*sigv
highb=5*sigv
+ write(iout,*) "lowb",lowb,highb,ii
d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+ if (itype(innt_org+i-1,4).gt.ntyp_molec(4)) d_t_work_new(ii)=d_t_work_new(ii)*1.0d-6
EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
+#ifdef DEBUG
write (iout,*) "i",i," ii",ii," geigen",geigen(i), &
" d_t_work_new",d_t_work_new(ii),innt_org+i-1
+#endif
enddo
enddo
- if (molnum(innt_org).ge.4) then
- mnum=molnum(innt_org)
+ if (molnum(innt_org+1).ge.4) then
+ mnum=molnum(innt_org+1)
do k=1,3
do i=1,n
ind=(i-1)*3+k
do j=1,n
d_t_work(ind)=d_t_work(ind)&
+Gvec(i,j)*d_t_work_new((j-1)*3+k)
+! if (Gvec(i,j).ne.0) write(iout,*) "for prot",ind,d_t_work(ind),d_t_work_new((j-1)*3+k),Gvec(i,j)
enddo
+#ifdef DEBUG
+ write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+#endif
enddo
enddo
endif
Im(2,2)=Im(2,2)+mp(mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
Im(3,3)=Im(3,3)+mp(mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
enddo
-
+ print *,"after pep",Im
do i=nnt,nct
mnum=molnum(i)
iti=iabs(itype(i,mnum))
Im(2,2)=Im(2,2)+msc(iabs(iti),mnum)*(pr(3)*pr(3)+pr(1)*pr(1))
Im(3,3)=Im(3,3)+msc(iabs(iti),mnum)*(pr(1)*pr(1)+pr(2)*pr(2))
enddo
+ print *,"after msc",Im
do i=nnt,nct-1
mnum=molnum(i)
mnum1=molnum(i+1)
Im(3,3)=Im(3,3)+Ip(mnum)*(1-dc_norm(3,i)*dc_norm(3,i))*&
vbld(i+1)*vbld(i+1)*0.25d0
enddo
+ print *,"after Ip",Im
do i=nnt,nct
mnum=molnum(i)
mnum1=molnum(i+1)
dc_norm(2,inres))*vbld(inres)*vbld(inres)
Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)*&
dc_norm(3,inres))*vbld(inres)*vbld(inres)
+ print *,i,inres,vbld(inres),iti,dc_norm(1,inres),Im(1,1)
endif
enddo
-
+ print *,"after ISC",Im
+ print *,"before angnom",L
+ print *,"before angnom",Im
call angmom(cm,L)
Im(2,1)=Im(1,2)
Im(3,1)=Im(1,3)
enddo
enddo
!c Finding the eigenvectors and eignvalues of the inertia tensor
+ print *,"before djacob",Imcp,eigvec,eigval
call djacob(3,3,10000,1.0d-10,Imcp,eigvec,eigval)
+ print *, "after djacob"
do i=1,3
if (dabs(eigval(i)).gt.1.0d-15) then
Id(i,i)=1.0d0/eigval(i)
M_PEP=M_PEP+mp(mnum)
do j=1,3
+ write(iout,*) "check",c(j,i),0.5d0*dc(j,i)
cm(j)=cm(j)+(c(j,i)+0.5d0*dc(j,i))*mp(mnum)
enddo
enddo
! do j=1,3
! cm(j)=mp(1)*cm(j)
! enddo
+ write(iout,*),"after pep",M_PEP,cm(1),cm(2),cm(3)
M_SC=0.0d0
do i=nnt,nct
mnum=molnum(i)
do j=1,3
cm(j)=cm(j)/(M_SC+M_PEP)
enddo
-! write(iout,*) "Center of mass:",cm
+ write(iout,*) "Center of mass:",cm
do i=nnt,nct-1
mnum=molnum(i)
! if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
endif
enddo
call angmom(cm,L)
-! write(iout,*) "The angular momentum after adjustment:"
-! write(iout,*) (L(j),j=1,3)
+ write(iout,*) "The angular momentum after adjustment:"
+ write(iout,*) (L(j),j=1,3)
return
end subroutine inertia_tensor
enddo
do i=nnt,nct-1
mnum=molnum(i)
- if (mnum.ge.5) mp(mnum)=msc(itype(i,mnum),mnum)
+ if (mnum.ge.4) mp(mnum)=msc(itype(i,mnum),mnum)
do j=1,3
pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
enddo
endif
! print *,i,pr,"pr",v
call vecpr(pr(1),v(1),vp)
-! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
-! " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
+ write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),&
+ " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
do j=1,3
L(j)=L(j)+mscab*vp(j)
enddo
! include 'COMMON.IOUNITS'
!el real(kind=8),dimension(6*nres) :: gamvec !(MAXRES6) maxres6=6*maxres
!el common /syfek/ gamvec
+!#define DEBUG
#ifdef FIVEDIAG
integer iposc,ichain,n,innt,inct
double precision rs(nres*2)
- real(kind=8) ::v_work(3,6*nres),vvec(2*nres)
+ real(kind=8) ::v_work(3,2*nres),vvec(2*nres)
#else
real(kind=8) :: v_work(6*nres)
#endif
nres6=6*nres
lprn=.false.
checkmode=.false.
+#ifdef FIVEDIAG
+ do i=1,nres2
+ do j=1,3
+ v_work(j,i)=0.0d0
+ enddo
+ enddo
+#else
+ do i=1,6*nres
+ v_work(i)=0.0d0
+ enddo
+#endif
! if (large) lprn=.true.
! if (large) checkmode=.true.
#ifdef FIVEDIAG
endif
#endif
return
+#undef DEBUG
end subroutine friction_force
!-----------------------------------------------------------------------------
subroutine setup_fricmat
inct=chain_border(2,ichain)
do i=innt,inct
mnum=molnum(i)
- if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+! if (iabs(itype(i,1)).ne.10.and.mnum.le.2) then
+ if (iabs(itype(i,1)).ne.10) then
ind=ind+2
else
DU1fric(ind)=gamvec(i-nnt+1)/4
#ifdef FIVEDIAG
integer ichain,innt,inct,iposc
#endif
-
+!#define DEBUG
do i=0,2*nres
do j=1,3
stochforc(j,i)=0.0d0
+ force(j,i)=0.0d0
enddo
enddo
x=0.0d0
#endif
! Compute the stochastic forces acting on bodies. Store in force.
do i=nnt,nct-1
+#ifdef FIVEDIAG
+ if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
+#endif
sig=stdforcp(i)
lowb=-5*sig
highb=5*sig
force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2)
enddo
enddo
+#ifdef DEBUG
+ write (iout,*) "Stochastic forces on sites"
+ do i=1,nres
+ write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3),&
+ (force(j,i+nres),j=1,3)
+ enddo
+#endif
+
#ifdef MPI
time_fsample=time_fsample+MPI_Wtime()-time00
#else
do j=1,3
stochforcvec(ind+j)=0.5d0*force(j,innt)
enddo
- if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt).ge.3) then
+ if (iabs(itype(innt,molnum(innt))).eq.10.or.molnum(innt+1).ge.3) then
do j=1,3
stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres)
enddo
endif
#endif
return
+#undef DEBUG
end subroutine stochastic_force
!-----------------------------------------------------------------------------
subroutine sdarea(gamvec)
allocate(DM(nres2),DU1(nres2),DU2(nres2))
allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2))
!#ifdef DEBUG
- allocate(Gvec(1300,1300))!maximum number of protein in one chain
+ allocate(Gvec(1300*2,1300*2))!maximum number of protein in one chain
!#endif
#else
write (iout,*) "Before A Allocation",nfgtasks-1
print *,"FIVEDIAG"
write (iout,*) "before FIVEDIAG"
#ifdef FIVEDIAG
- if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
+ if(.not.allocated(Ghalf)) allocate(Ghalf(1300*(1300+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
#endif
#ifndef FIVEDIAG
if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2)
real(kind=8),dimension(:),allocatable :: xsolv,rs
real(kind=8),dimension(:,:),allocatable :: accel
integer i,j,ind,ichain,n,iposc,innt,inct,inct_prev,mnum
- accel=0.0d0
+! accel=0.0d0
!#define DEBUG
if (.not.allocated(forces)) allocate(forces(6*nres))
if (.not.allocated(d_a_vec)) allocate(d_a_vec(6*nres))
- if (.not.allocated(xsolv)) allocate(xsolv(3*ndim))
- if (.not.allocated(rs)) allocate(rs(3*ndim))
+ if (.not.allocated(xsolv)) allocate(xsolv(ndim))
+ if (.not.allocated(rs)) allocate(rs(ndim))
if (.not.allocated(accel)) allocate(accel(3,0:2*nres))
+ accel=0.0d0
#ifdef DEBUG
do i=1,6*nres
ind=(i-1)*3+k+1
d_a_tmp(ind)=0.0d0
do j=1,dimen
-! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1
+! write (2,*) "k,i,j,ind",k,i,j,ind,(j-1)*3+k+1,&
! call flush(2)
-! & Ginv(i,j),z((j-1)*3+k+1),
-! & Ginv(i,j)*z((j-1)*3+k+1)
+! Ginv(i,j),z((j-1)*3+k+1), &
+! Ginv(i,j)*z((j-1)*3+k+1)
d_a_tmp(ind)=d_a_tmp(ind) &
+Ginv(j,i)*z((j-1)*3+k+1)
! d_a_tmp(ind)=d_a_tmp(ind)
endif
ncont=0
kkk=0
-! print *,'nnt=',nnt,' nct=',nct
+ print *,'nnt=',nnt,' nct=',nct,nnt_molec(1),nct_molec(1)-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))
xmedi=xi+0.5*dxi
ymedi=yi+0.5*dyi
zmedi=zi+0.5*dzi
- do 4 j=i+2,nct-1
+ do 4 j=i+2,nct_molec(1)-1
if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) goto 4
iteli=itel(i)
itelj=itel(j)
do i=nz_start,nz_end
! iatom=iatom+1
iti=itype(i+nstart_seq-nstart_sup,molnum(i))
- if (iti.eq.ntyp1_molec(molnum(i))) cycle
+ if (iti.gt.ntyp_molec(molnum(i))) cycle
+ if ((molnum(i).eq.5).and.(itype(i,5).eq.-1)) cycle
iatom=iatom+1
mnum=molnum(i+nstart_seq-nstart_sup)
do k=1,3
nct_molec(i)=itmp
endif
enddo
+ if (itype(1,molnum(1)).eq.ntyp1_molec(1)) then
+ nnt_molec(1)=2
+ else
+ nnt_molec(1)=1
+ endif
! nct_molec(1)=nres_molec(1)-1
itmp=0
do i=2,5
-sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
!el allocate(jgrad_start(igrad_start:igrad_end))
!el allocate(jgrad_end(igrad_start:igrad_end)) !(maxres)
- jgrad_start(igrad_start)= &
- ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2 &
- +igrad_start
- jgrad_end(igrad_start)=nres
- if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
- jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2 &
- +igrad_end
- do i=igrad_start+1,igrad_end-1
- jgrad_start(i)=i+1
- jgrad_end(i)=nres
- enddo
+! jgrad_start(igrad_start)= &
+! ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2 &
+! +igrad_start
+! jgrad_end(igrad_start)=nres
+! if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
+! jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2 &
+! +igrad_end
+! do i=igrad_start+1,igrad_end-1
+! jgrad_start(i)=i+1
+! jgrad_end(i)=nres
+! enddo
! THIS SHOULD BE FIXED
itmp=0
do i=1,4
! include 'COMMON.SETUP'
#ifdef MPI
call int_bounds(nhpb,link_start,link_end)
+ call int_bounds(cnhpb,clink_start,clink_end)
write (iout,*) 'Processor',fg_rank,' CG group',kolor,&
' absolute rank',MyRank,&
' nhpb',nhpb,' link_start=',link_start,&
#else
link_start=1
link_end=nhpb
+ clink_start=1
+ clink_end=cnhpb
+
#endif
return
end subroutine hpb_partition
!-----------------------------------------------------------------------------
logical preminim, forceminim ! pre-minimizaation flag
#ifdef FIVEDIAG
- integer,parameter :: maxchain=50
+ integer,parameter :: maxchain=1000
integer,dimension(maxchain) :: dimen_chain,iposd_chain
real(kind=8),dimension(:),allocatable :: DMfric,DU1fric,&
DU2fric
! common /cntrl/
integer :: modecalc,iscode,indpdb,indback,indphi,iranconf,&
icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode,&
- tubemode,genconstr,afmlog,selfguide,scelemode,tor_mode,oldion,DRY_MARTINI
+ tubemode,genconstr,afmlog,selfguide,scelemode,tor_mode,oldion,DRY_MARTINI,&
+ vacuum
logical :: minim,refstr,pdbref,overlapsc,&
energy_dec,sideadd,lsecondary,read_cart,unres_pdb,&
vdisulf,searchsc,lmuca,dccart,extconf,out1file,&
integer,dimension(:),allocatable :: iss !(maxss)
! common /links/
real(kind=8),dimension(:),allocatable :: dhpb,forcon,dhpb1,fordepth !(maxdim) !el dhpb1 !!! nie używane
- integer :: nhpb
- integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
+ real(kind=8),dimension(:),allocatable :: cdhpb,cforcon,cdhpb1,cfordepth !(maxdim) !el dhpb1 !!! nie używane
+
+ integer :: nhpb,cnhpb
+ integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb,cihpb,cjhpb,cibecarb !(maxdim) !el ibecarb !!! nie używane
! common /restraints/
real(kind=8) :: weidis
! common /links_split/
- integer :: link_start,link_end
+ integer :: link_start,link_end,clink_start,clink_end
! common /dyn_ssbond/
real(kind=8) :: Ht,atriss,btriss,ctriss,dtriss
integer,dimension(:),allocatable :: idssb,jdssb !(maxdim)
call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
endif
if (nres_molec(1).gt.0) then
+! print *,"before make_SCp_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
-! write (iout,*) "after make_SCp_inter_list"
+! print *, "after make_SCp_inter_list"
if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
! write (iout,*) "after make_SCSC_inter_list"
if (nres_molec(4).gt.0) then
if (mod(itime_mat,imatupdate).eq.0) then
! print *,'Processor',myrank,' calling etotal ipot=',ipot
call make_cat_pep_list
+! write(iout,*) "after make cat pep"
! call make_cat_cat_list
endif
endif
! print *,'Processor',myrank,' calling etotal ipot=',ipot
! call make_cat_pep_list
call make_cat_cat_list
+ write(iout,*) "after make cat cat"
endif
endif
! write (iout,*) "after make_pp_inter_list"
ehpb=0.0D0
! write(iout,*)'edis: nhpb=',nhpb!,' fbr=',fbr
! write(iout,*)'link_start=',link_start,' link_end=',link_end
- if (link_end.eq.0) return
- do i=link_start,link_end
+ if (clink_end.eq.0) return
+ do i=clink_start,clink_end
! If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
! CA-CA distance used in regularization of structure.
- ii=ihpb(i)
- jj=jhpb(i)
+ ii=cihpb(i)
+ jj=cjhpb(i)
! iii and jjj point to the residues for which the distance is assigned.
if (ii.gt.nres) then
iii=ii-nres
dd=dist(ii,jj)
if (constr_dist.eq.11) then
ehpb=ehpb+fordepth(i)**4.0d0 &
- *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ *rlornmr1(dd,cdhpb(i),cdhpb1(i),cforcon(i))
fac=fordepth(i)**4.0d0 &
- *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ *rlornmr1prim(dd,cdhpb(i),cdhpb1(i),cforcon(i))/dd
if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
- ehpb,fordepth(i),dd
+ ehpb,cfordepth(i),dd
else
- if (dhpb1(i).gt.0.0d0) then
- ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
- fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+ if (cdhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*cforcon(i)*gnmr1(dd,cdhpb(i),cdhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
!c write (iout,*) "beta nmr",
!c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
else
dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ rdis=dd-cdhpb(i)
!C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=cforcon(i)
!C Calculate the contribution to energy.
ehpb=ehpb+waga*rdis*rdis
!c write (iout,*) "beta reg",dd,waga*rdis*rdis
dd=dist(ii,jj)
if (constr_dist.eq.11) then
- ehpb=ehpb+fordepth(i)**4.0d0 &
- *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
+ ehpb=ehpb+cfordepth(i)**4.0d0 &
+ *rlornmr1(dd,dhpb(i),cdhpb1(i),cforcon(i))
fac=fordepth(i)**4.0d0 &
- *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
+ *rlornmr1prim(dd,cdhpb(i),cdhpb1(i),cforcon(i))/dd
if (energy_dec) write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, &
ehpb,fordepth(i),dd
else
- if (dhpb1(i).gt.0.0d0) then
- ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
- fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+ if (cdhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*cforcon(i)*gnmr1(dd,cdhpb(i),cdhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,cdhpb(i),cdhpb1(i))/dd
!c write (iout,*) "alph nmr",
!c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
else
vec(2)=yj
vec(3)=zj
dd=sqrt(xj*xj+yj*yj+zj*zj)
- rdis=dd-dhpb(i)
+ rdis=dd-cdhpb(i)
!C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=cforcon(i)
!C Calculate the contribution to energy.
ehpb=ehpb+waga*rdis*rdis
- if (energy_dec) write (iout,'(a6,2i5,5f10.3)') "edis",ii,jj, &
- ehpb,dd,dhpb(i),waga,rdis
+ ! if (energy_dec)
+! write (iout,'(a6,2i5,5f10.3)') "edisip",ii,jj, &
+! ehpb,dd,cdhpb(i),waga,rdis
!c write (iout,*) "alpha reg",dd,waga*rdis*rdis
!C
#ifdef MPI
include 'mpif.h'
#endif
- real(kind=8),dimension(3,-1:nres) :: gradbufc,gradbufx,gradbufc_sum,&
+ real(kind=8),dimension (:,:), allocatable :: gradbufc,gradbufx,gradbufc_sum,&
gloc_scbuf !(3,maxres)
- real(kind=8),dimension(4*nres) :: glocbuf !(4*maxres)
+ real(kind=8),dimension(:),allocatable :: glocbuf !(4*maxres)
!#endif
!el local variables
integer :: i,j,k,ierror,ierr
gscloc_norm,gvdwx_norm,gradx_scp_norm,ghpbx_norm,&
gradxorr_norm,gsccorrx_norm,gsclocx_norm,gcorr6_max,&
gsccorr_max,gsccorrx_max,time00
+ if (.not.allocated(gradbufc)) then
+ allocate(gradbufc(3,-1:nres))
+ allocate(gradbufx(3,-1:nres))
+ allocate(gradbufc_sum(3,-1:nres))
+ allocate(gloc_scbuf(3,-1:nres))
+ allocate(glocbuf(4*nres))
+ endif
+! print *,"in sum gradient"
! include 'COMMON.SETUP'
! include 'COMMON.IOUNITS'
! include 'COMMON.FFIELD'
! enddo
! enddo
!el#define DEBUG
+! print *,"gradbufc after summing"
#ifdef DEBUG
write (iout,*) "gradbufc after summing"
do i=1,nres
!el ind=ind+1
itypj=iabs(itype(j,1))
if (itypj.eq.ntyp1) cycle
+ if (vacuum.gt.0) then
+ if (itypi.eq.10) cycle
+ if (itypj.eq.10) cycle
+ endif
+
CALL elgrad_init(eheadtail,Egb,Ecl,Elj,Equad,Epol)
! if (j.ne.78) cycle
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) :: dist_init, dist_temp,r_buff_list
- integer:: contlisti(250*nres),contlistj(250*nres)
+ integer:: contlisti(250*nres_molec(1)),contlistj(250*nres_molec(1))
! integer :: newcontlisti(200*nres),newcontlistj(200*nres)
integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_sc,g_ilist_sc
integer displ(0:nprocs),i_ilist_sc(0:nprocs),ierr
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) :: dist_init, dist_temp,r_buff_list
- integer:: contlistscpi(350*nres),contlistscpj(350*nres)
+ integer:: contlistscpi(350*nres_molec(1)),contlistscpj(350*nres_molec(1))
! integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr
! print *,"START make_SC"
r_buff_list=5.0
ilist_scp=0
+ print *,"in make SCp list", iatscp_s,iatscp_e
do i=iatscp_s,iatscp_e
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
xi=0.5D0*(c(1,i)+c(1,i+1))
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
- integer:: contlistppi(250*nres),contlistppj(250*nres)
+ integer:: contlistppi(250*nres_molec(1)),contlistppj(250*nres_molec(1))
! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
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(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)
- integer:: contlistcatscangi(250*nres),contlistcatscangj(250*nres)
- integer:: contlistcatscangfi(250*nres),contlistcatscangfj(250*nres),&
- contlistcatscangfk(250*nres)
- integer:: contlistcatscangti(250*nres),contlistcatscangtj(250*nres)
- integer:: contlistcatscangtk(250*nres),contlistcatscangtl(250*nres)
+ integer,dimension(:), allocatable :: contlistcatpnormi,&
+ contlistcatpnormj, contlistcatscnormi,contlistcatscnormj,&
+ contlistcatptrani,contlistcatptranj,contlistcatsctrani,&
+ contlistcatsctranj,contlistcatscangi,contlistcatscangj,&
+ contlistcatscangfi,contlistcatscangfj,contlistcatscangfk,&
+ contlistcatscangti,contlistcatscangtj,contlistcatscangtk,&
+ contlistcatscangtl
! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
ilist_catptran=0
ilist_catsctran=0
ilist_catscang=0
+ if (.not.allocated(contlistcatpnormi)) then
+ allocate(contlistcatpnormi(300*nres),contlistcatpnormj(300*nres))
+ allocate(contlistcatscnormi(250*nres),contlistcatscnormj(250*nres))
+ allocate(contlistcatptrani(250*nres),contlistcatptranj(250*nres))
+ allocate(contlistcatsctrani(250*nres),contlistcatsctranj(250*nres))
+ allocate(contlistcatscangi(250*nres),contlistcatscangj(250*nres))
+ allocate(contlistcatscangfi(250*nres),contlistcatscangfj(250*nres),&
+ contlistcatscangfk(250*nres))
+ allocate(contlistcatscangti(250*nres),contlistcatscangtj(250*nres))
+ allocate(contlistcatscangtk(250*nres),contlistcatscangtl(250*nres))
+ endif
r_buff_list=6.0
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:: contlistmartpi(300*nres),contlistmartpj(300*nres)
- integer:: contlistmartsci(250*nres),contlistmartscj(250*nres)
+ integer,dimension(:),allocatable:: contlistmartpi,contlistmartpj,&
+ contlistmartsci,contlistmartscj
! integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
! write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
ilist_martp=0
ilist_martsc=0
-
+ if (.not.allocated(contlistmartpi)) then
+ allocate(contlistmartpi(300*nres),contlistmartpj(300*nres))
+ allocate(contlistmartsci(250*nres),contlistmartscj(250*nres))
+ endif
r_buff_list=6.0
itmp=0
double precision energia(0:n_ene)
double precision x(nvar),g(nvar)
integer i
+ print *,"in funcgrad"
+! if (not.allocated(x)) then
+! allocate x(nvar)
+! allocate g(nvar)
+! endif
call var_to_geom(nvar,x)
call zerograd
call chainbuild
call etotal(energia(0))
+ print *,"before sum gradient"
call sum_gradient
funcgrad=energia(0)
+ print *,"before cart2intgrad"
call cart2intgrad(nvar,g)
if (usampl) then
do i=1,nres-3
do i=1,nvar
g(i)=g(i)+gloc(i,icg)
enddo
+ print *,"finish funcgrad"
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,dimension(:,:,:),allocatable :: drt,rdt,prordt,prodrt
+ double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp,dp(3,3),temp(3,3)
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
+ print *,"in cart2intgrad"
+ if (.not.allocated(drt)) then
+ allocate(drt(3,3,nres))
+ allocate(rdt(3,3,nres))
+ allocate(prordt(3,3,nres))
+ allocate(prodrt(3,3,nres))
+ endif
sideonly=.false.
g=0.0d0
if (sideonly) goto 10
! dyj = dc_norm( 2, nres+j )
! dzj = dc_norm( 3, nres+j )
- itypi = itype(i,1)
+ itypi = iabs(itype(i,1))
itypj = itype(j,4)
! Parameters from fitting the analitical expressions to the PMF obtained by umbrella
! sampling performed with amber package
IF (rij_shift.le.0.0D0) THEN
evdw = 1.0D20
if (evdw.gt.1.0d6) then
- write (*,'(2(1x,a3,i3),7f7.2)') &
- restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
- 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq
- write(*,*) facsig,faceps1_inv,om1,chiom1,chi1
- write(*,*) "ANISO?!",chi1
+! write (*,'(2(1x,a3,i3),7f7.2)') &
+! restyp(itype(i,1),1),i,restyp(itype(j,1),1),j,&
+! 1.0d0/rij,Rtail,Rhead,rij_shift, sig, sig0ij,sigsq
+! write(*,*) facsig,faceps1_inv,om1,chiom1,chi1
+! write(*,*) "ANISO?!",chi1
!evdwij,Fcav,Ecl,Egb,Epol,Fisocav,Elj,&
! Equad,evdwij+Fcav+eheadtail,evdw
endif
facd2 = dtailmart(2,itypi,itypj) * vbld_inv(j)
DO k = 1, 3
pom = ertail(k)-facd1*(ertail(k)-erdxi*dC_norm(k,i+nres))
+! if ((i.eq.40).and.(j.eq.376)) write(iout,*) "before grad",gradpepmartx(k,i)
+! if ((i.eq.40).and.(j.eq.376)) write(iout,*) "during grad", dFdR ,gg(k),sss_ele_cut,(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
gradpepmartx(k,i) = gradpepmartx(k,i) &
- (( dFdR + gg(k) ) * pom)*sss_ele_cut&
-(evdwij+Fcav)*rij*sss_ele_grad*rreal(k)
ENDDO
!c! Compute head-head and head-tail energies for each state
!! if (.false.) then ! turn off electrostatic
+! if ((i.eq.40).and.(j.eq.376)) write(iout,*) "after grad evdw",(gradpepmartx(k,i),k=1,3)
+
isel = iabs(Qi)+iabs(Qj)
if ((itype(j,4).gt.4).and.(itype(j,4).lt.14)) isel=isel+2
! isel=0
ENDIF
! write(iout,*) "not yet implemented",j,itype(j,5)
!! endif ! turn off electrostatic
+! if ((i.eq.40).and.(j.eq.376)) write(iout,*) "after grad elec",(gradpepmartx(k,i),k=1,3)
evdw = evdw + (Fcav + eheadtail)*sss_ele_cut
! if (evdw.gt.1.0d6) then
! write (*,'(2(1x,a3,i3),3f6.2,10f16.7)') &
! print *,"before sc_grad_mart", i,j, gradpepmart(1,j)
! iF (nstate(itypi,itypj).eq.1) THEN
CALL sc_grad_mart
+! if ((i.eq.40).and.(j.eq.376)) write(iout,*) "after sc_grad ",(gradpepmartx(k,i),k=1,3)
! print *,"after sc_grad_mart", i,j, gradpepmart(1,j)
! END IF
use calc_data
real(kind=8) :: eheadtail,Egb,Ecl,Elj,Equad,Epol,Rb
eps_out=80.0d0
- itypi = itype(i,1)
+ itypi = iabs(itype(i,1))
itypj = itype(j,4)
! print *,"in elegrad",i,j,itypi,itypj
!c! 1/(Gas Constant * Thermostate temperature) = BetaT
return
end subroutine alloc_geo_arrays
!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
subroutine returnbox
+ real(kind=8),dimension(3) :: boxx,vecc
+ integer :: i,j,k
+ boxx(1)=boxxsize
+ boxx(2)=boxysize
+ boxx(3)=boxzsize
+ do i=1,nres
+ if (i.gt.2) then
+ if (molnum(i-1).ne.molnum(i)) then
+ do j=1,3
+ vecc(j)=c(j,i)-dmod(c(j,i),boxx(j))-(boxshift22(dmod(c(j,i),boxx(j))-c(j,2),boxx(j)))
+ enddo
+ endif
+ endif
+ if (molnum(i).gt.2) then
+ do j=1,3
+ c(j,i)=c(j,i)-vecc(j)
+! c(j,i)=c(j,i)+(boxshift22(c(j,i)-c(j,2),boxx(j)))
+! if (dabs(c(j,i)-c(j,2)).gt.boxx(j)/2) print *,"error",i,j,c(j,i),c(j,2),vecc(j)
+! if (c(j,i).lt.0) c(j,i)=c(j,i)+boxx(j)
+! if (c(j,2).lt.0) c(j,i)=c(j,i)-boxx(j)
+ enddo
+ else
+ if (i.eq.2) then
+ do j=1,3
+ vecc(j)=c(j,i)-dmod(c(j,i),boxx(j))
+ enddo
+ endif
+ do j=1,3
+ c(j,i)=c(j,i)-vecc(j)
+ c(j,i+nres)=c(j,i+nres)-vecc(j)
+ enddo
+ endif
+ enddo
+ return
+ end subroutine returnbox
+
+! if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
+! print *,"return box,before wrap",c(1,i)
+
+!-----------------------------------------------------------------------------
+ subroutine returnbox2
integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
integer :: chain_end,ireturnval,idum,mnumi1
real*8 :: difference,xi,boxsize,x,xtemp,box2shift
k=k+1
endif
endif
-! print *,"tu2",i,k
+ print *,"tu2",i,k
if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1
if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1
-! print *,"tu2",i,k
+ print *,"tu2",i,k
do j=1,3
c(j,i)=dmod(c(j,i),boxx(j))
box2shift=xtemp
endif
xi=box2shift
-! print *,c(1,2),xi,xtemp,box2shift
+ print *,c(1,2),xi,xtemp,box2shift
c(j,i)=c(j,k)+xi
enddo
do j=1,3
! enddo
! endif
return
- end subroutine returnbox
+ end subroutine returnbox2
+ double precision function boxshift22(x,boxsize)
+ implicit none
+ double precision x,boxsize
+ double precision xtemp
+ xtemp=dmod(x,boxsize)
+ if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+ boxshift22=-boxsize
+ else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+ boxshift22=boxsize
+ else
+ boxshift22=0
+ endif
+ return
+ end function boxshift22
+
!-------------------------------------------------------------------------------------------------------
end module geometry
chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
enddo
allocate(tabpermchain(50,5040))
+ if (symetr.eq.1) then
+ npermchain=1
+ tabpermchain(1,1)=1
+ else
call chain_symmetry(npermchain,tabpermchain)
+ endif
print *,'NNT=',NNT,' NCT=',NCT
if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
print *,"in seq2"
ichain=1
new_chain=.true.
- if (.not.allocated(chain_length)) allocate(chain_length(50))
- if (.not.allocated(chain_border)) allocate(chain_border(2,50))
+ if (.not.allocated(chain_length)) allocate(chain_length(1000))
+ if (.not.allocated(chain_border)) allocate(chain_border(2,1000))
chain_length(ichain)=0
ii=1
integer itemp(50),&
npermchain,tabpermchain(50,5040),&
tabperm(50,5040),mapchain(50),&
- iequiv(50,nres),iflag(nres)
+ iflag(nres)
integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
nperm,npermc,ind,mnum
+ integer,dimension(:,:),allocatable :: iequiv
+ if (.not.allocated(iequiv)) allocate(iequiv(1000,nres))
if (nchain.eq.1) then
npermchain=1
tabpermchain(1,1)=1
! write (iout,*) "Calling read_dist_constr"
! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
! call flush(iout)
- if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
- if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
- if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
- if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
- if(.not.allocated(forcon)) allocate(forcon(maxdim))
- if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
- if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim))
+ if(.not.allocated(cihpb)) allocate(cihpb(maxdim))
+ if(.not.allocated(cjhpb)) allocate(cjhpb(maxdim))
+ if(.not.allocated(cdhpb)) allocate(cdhpb(maxdim))
+ if(.not.allocated(cdhpb1)) allocate(cdhpb1(maxdim))
+ if(.not.allocated(cforcon)) allocate(cforcon(maxdim))
+ if(.not.allocated(cfordepth)) allocate(cfordepth(maxdim))
+ if(.not.allocated(cibecarb)) allocate(cibecarb(maxdim))
if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
call gen_dist_constr2
go to 1712
endif
+ cnhpb=0
call card_concat(controlcard,.true.)
call readi(controlcard,"NFRAG",nfrag_,0)
call readi(controlcard,"NPAIR",npair_,0)
write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
ddjk=dist(j,k)
if (constr_dist.eq.1) then
- nhpb=nhpb+1
- ihpb(nhpb)=j
- jhpb(nhpb)=k
- dhpb(nhpb)=ddjk
- forcon(nhpb)=wfrag_(i)
+ cnhpb=cnhpb+1
+ cihpb(cnhpb)=j
+ cjhpb(cnhpb)=k
+ cdhpb(cnhpb)=ddjk
+ cforcon(cnhpb)=wfrag_(i)
else if (constr_dist.eq.2) then
if (ddjk.le.dist_cut) then
- nhpb=nhpb+1
- ihpb(nhpb)=j
- jhpb(nhpb)=k
- dhpb(nhpb)=ddjk
- forcon(nhpb)=wfrag_(i)
+ print *,"tu",nhpb,i,j,k,maxdim
+ cnhpb=cnhpb+1
+ cihpb(cnhpb)=j
+ cjhpb(cnhpb)=k
+ cdhpb(cnhpb)=ddjk
+ cforcon(cnhpb)=wfrag_(i)
endif
else
- nhpb=nhpb+1
- ihpb(nhpb)=j
- jhpb(nhpb)=k
- dhpb(nhpb)=ddjk
- forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ cnhpb=cnhpb+1
+ cihpb(cnhpb)=j
+ cjhpb(cnhpb)=k
+ cdhpb(cnhpb)=ddjk
+ cforcon(cnhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king) &
write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
#else
write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
#endif
enddo
enddo
endif
do j=ifrag_(1,ii),ifrag_(2,ii)
do k=ifrag_(1,jj),ifrag_(2,jj)
- nhpb=nhpb+1
- ihpb(nhpb)=j
- jhpb(nhpb)=k
- forcon(nhpb)=wpair_(i)
- dhpb(nhpb)=dist(j,k)
+ cnhpb=cnhpb+1
+ cihpb(cnhpb)=j
+ cjhpb(cnhpb)=k
+ cforcon(cnhpb)=wpair_(i)
+ cdhpb(cnhpb)=dist(j,k)
#ifdef MPI
if (.not.out1file .or. me.eq.king) &
write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ nhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
#else
write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ nhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
#endif
enddo
enddo
! nhpb=nhpb+1
! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
if (constr_dist.eq.11) then
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
- ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
+ read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
+ cibecarb(i),cforcon(cnhpb+1),cfordepth(cnhpb+1)
!EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
- fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0)
- forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0)
+ fordepth(nhpb+1)=cfordepth(cnhpb+1)**(0.25d0)
+ forcon(nhpb+1)=cforcon(cnhpb+1)**(0.25d0)
else
!C print *,"in else"
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
- ibecarb(i),forcon(nhpb+1)
+ read (inp,*) cihpb(cnhpb+1),cjhpb(cnhpb+1),cdhpb(i),cdhpb1(i), &
+ cibecarb(i),cforcon(cnhpb+1)
endif
- if (forcon(nhpb+1).gt.0.0d0) then
- nhpb=nhpb+1
+ if (cforcon(nhpb+1).gt.0.0d0) then
+ cnhpb=cnhpb+1
if (ibecarb(i).gt.0) then
- ihpb(i)=ihpb(i)+nres
- jhpb(i)=jhpb(i)+nres
+ cihpb(i)=ihpb(i)+nres
+ cjhpb(i)=jhpb(i)+nres
endif
- if (dhpb(nhpb).eq.0.0d0) &
- dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ if (cdhpb(nhpb).eq.0.0d0) &
+ cdhpb(cnhpb)=dist(ihpb(cnhpb),jhpb(cnhpb))
endif
#ifdef MPI
if (.not.out1file .or. me.eq.king) &
write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
#else
write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(cnhpb),cjhpb(cnhpb),cdhpb(cnhpb),cforcon(cnhpb)
#endif
enddo
1712 continue
do j=i+2,nstart_sup+nsup-1
distance=dist(i,j)
if (distance.le.15.0) then
- nhpb=nhpb+1
- ihpb(nhpb)=i+nstart_seq-nstart_sup
- jhpb(nhpb)=j+nstart_seq-nstart_sup
- forcon(nhpb)=sqrt(0.04*distance)
- fordepth(nhpb)=sqrt(40.0/distance)
- dhpb(nhpb)=distance-0.1d0
- dhpb1(nhpb)=distance+0.1d0
+ cnhpb=cnhpb+1
+ cihpb(nhpb)=i+nstart_seq-nstart_sup
+ cjhpb(nhpb)=j+nstart_seq-nstart_sup
+ cforcon(nhpb)=sqrt(0.04*distance)
+ cfordepth(nhpb)=sqrt(40.0/distance)
+ cdhpb(nhpb)=distance-0.1d0
+ cdhpb1(nhpb)=distance+0.1d0
#ifdef MPI
if (.not.out1file .or. me.eq.king) &
write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
#else
write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
- nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ cnhpb,cihpb(nhpb),cjhpb(nhpb),cdhpb(nhpb),cforcon(nhpb)
#endif
endif
enddo
else
do i=nstart_sup,nstart_sup+nsup-1
do j=i+2,nstart_sup+nsup-1
- nhpb=nhpb+1
- ihpb(nhpb)=i+nstart_seq-nstart_sup
- jhpb(nhpb)=j+nstart_seq-nstart_sup
- forcon(nhpb)=weidis
- dhpb(nhpb)=dist(i,j)
+ cnhpb=cnhpb+1
+ cihpb(nhpb)=i+nstart_seq-nstart_sup
+ cjhpb(nhpb)=j+nstart_seq-nstart_sup
+ cforcon(nhpb)=weidis
+ cdhpb(nhpb)=dist(i,j)
enddo
enddo
endif
if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
do i=1,ntyp
+ if (vacuum.gt.0) then
+ if (i.eq.10) cycle
+ endif
do j=1,i
-! write (*,*) "Im in ALAB", i, " ", j
+ if (vacuum.gt.0) then
+ if (j.eq.10) cycle
+ endif
+ write (*,*) "Im in ALAB", i, " ", j
read(isidep,*) &
eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
(alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
! enddo
! endif
! endif
-
+ do i=1,5
+ nres_molec(i)=0
+ enddo
+ do i=1,nres
+ nres_molec(molnum(i))=nres_molec(molnum(i))+1
+ enddo
+ print *,"MY FINAL NRES",nres,nres_molec
if (lprn) then
write (iout,'(/a)') &
"Cartesian coordinates of the reference structure after sorting"
(c(j,ires+nres),j=1,3)
enddo
endif
-
+
print *,seqalingbegin,nres
if(.not.allocated(vbld)) then
allocate(vbld(2*nres))
write (iout,*) "with_theta_constr ",with_theta_constr
AFMlog=(index(controlcard,'AFM'))
selfguide=(index(controlcard,'SELFGUIDE'))
+ vacuum=(index(controlcard,'VACUUM'))
print *,'AFMlog',AFMlog,selfguide,"KUPA"
call readi(controlcard,'GENCONSTR',genconstr,0)
call reada(controlcard,'BOXX',boxxsize,100.0d0)
COS45 = .70710678
S45SQ = 0.50
C45SQ = 0.50
+ print *,"in math",N,NMAX
! UNIT EIGENVECTOR MATRIX
DO 70 I = 1,N
DO 7 J = I,N
1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',&
'MENT = ',1PE14.7)
! ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER
+ print *,"before first AJJ?"
300 DO 14 I=1,N
14 AJJ(I)=A(I,I)
LT=N+1
17 AIIMIN=AJJ(I)
IT=I
16 CONTINUE
+ print *,IT,"IT"
IN=L
AII(IN)=AIIMIN
AJJ(IT)=1.0E+30
real(kind=8),DIMENSION(N,8) :: WRK
real(kind=8),DIMENSION(N) :: EIG
integer,DIMENSION(N) :: IWRK
- real(kind=8),DIMENSION(LDVECT,NVECT) :: VECTOR
+ real(kind=8),DIMENSION(:,:),allocatable :: VECTOR
!
!el integer :: IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
!el integer :: KDIAG,ICORFL,IXDR
! IERR = ERROR FLAG (OUTPUT)
! IWRK = N INTEGER WORDS OF SCRATCH SPACE
!
+ if (.not.allocated(VECTOR)) allocate(VECTOR(LDVECT,NVECT))
IERR = 0
!
! ----- USE STEVE ELBERT'S ROUTINE -----
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(:),allocatable :: x,d,xx !(maxvar) (maxvar=6*maxres)
+ real(kind=8),dimension(:),allocatable :: x,d,xx,g !(maxvar) (maxvar=6*maxres)
!el common /przechowalnia/ v
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
+ real(kind=8) :: f1,etot,rdum(1) !,fdum
#ifdef LBFGS
external optsave
maxiter=maxmin
v(32)=rtolf
! controls initial step size
v(35)=1.0D-1
+ print *,"before d components"
! large vals of d correspond to small components of step
+ print *,"before x allocate"
+ if (.not.allocated(x)) then
+ allocate(x(6*nres))
+ allocate(xx(6*nres))
+ allocate(d(6*nres))
+ allocate(g(6*nres))
+ endif
do i=1,6*nres
d(i)=1.0D-1
enddo
#endif
+ print *,"before x allocate"
if (.not.allocated(x)) then
allocate(x(6*nres))
allocate(xx(6*nres))
allocate(d(6*nres))
+ allocate(g(6*nres))
endif
+! print *,"after x allocate"
k=0
do i=1,nres-1
do j=1,3
enddo
enddo
do i=2,nres-1
- if (ialph(i,1).gt.0) then
+
+ if ((ialph(i,1).gt.0).and.(molnum(i).le.3)) then
do j=1,3
k=k+1
x(k)=dc(j,i+nres)
enddo
nvarx=k
call chainbuild_cart
+! print *,"before etotal"
call etotal(energia(0))
call enerprint(energia(0))
#ifdef LBFGS
end do
write (iout,*) "minim_dc Calling lbfgs"
call lbfgs (nvarx,x,etot,grdmin,funcgrad_dc,optsave)
+ print *,"after lbfgs"
deallocate(scale)
#else
! call check_ecartint
include 'mpif.h'
#endif
integer k,nf,i,j
- real(kind=8),dimension(nres*6):: x,g
+ real(kind=8),dimension(6*nres):: x,g
double precision energia(0:n_ene)
! common /gacia/ nf
!c
+! if (.not.allocated(x)) allocate(x(6*nres),g(6*nres))
nf=nf+1
k=0
do i=1,nres-1
call chainbuild_cart
call zerograd
call etotal(energia(0))
-!c write (iout,*) "energia",energia(0)
+! write (iout,*) "energia",energia(0)
funcgrad_dc=energia(0)
call cartgrad
k=0
!c
!c initialize some values to be used below
!c
+ print *,"in lbfgs",nvar
ncalls = 0
rms = sqrt(dble(nvar))
if (coordtype .eq. 'CARTESIAN') then
if (jprint .lt. 0) jprint = 1
if (iwrite .lt. 0) iwrite = 1
write(iout,*) "maxiter",maxiter
+ jprint=1
!c
!c set default parameters for the line search
!c
niter = nextiter - 1
maxiter = niter + maxiter
ncalls = ncalls + 1
+ print *,"just before fgvalue"
! itime_mat=niter
f = fgvalue (x0,g)
+ print *,"just after fgvalue"
f_old = f
m = 0
gamma = 1.0d0
!
! get new function and projected gradient following a step
!
+ print *,"my calls",ncalls
ncalls = ncalls + 1
! 3/14/2020 Adam Liwo: added the condition to prevent from infinite
! iteration loop
allocate(sequenc(maxres))
allocate(molnum(maxres))
allocate(itype(maxres,5))
+ allocate(ihpb(100))
+ allocate(jhpb(100))
+
if (iargc().lt.6) then
print '(2a)',&
"Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",&
projects / unres4.git / commitdiff