From: Adam Sieradzan Date: Thu, 18 Jan 2024 14:02:30 +0000 (+0100) Subject: pleanty of changes X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=f3be969fc76e5a94d38a3dbc5cadd8c95a27067d;p=unres4.git pleanty of changes --- diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt index d9be5ea..41d35fa 100644 --- a/source/unres/CMakeLists.txt +++ b/source/unres/CMakeLists.txt @@ -82,7 +82,7 @@ if (Fortran_COMPILER_NAME STREQUAL "ifort") 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 ") diff --git a/source/unres/MD.F90 b/source/unres/MD.F90 index bf993eb..52f96af 100644 --- a/source/unres/MD.F90 +++ b/source/unres/MD.F90 @@ -789,13 +789,13 @@ 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) @@ -808,7 +808,7 @@ 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) @@ -824,7 +824,7 @@ 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. @@ -850,7 +850,7 @@ 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) @@ -859,7 +859,7 @@ 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 @@ -869,12 +869,12 @@ 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 @@ -886,16 +886,16 @@ 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 !----------------------------------------------------------------------------- @@ -1002,8 +1002,8 @@ 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 @@ -1255,7 +1255,7 @@ 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 @@ -1538,6 +1538,8 @@ ! 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)) & @@ -2290,10 +2292,12 @@ ! 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 @@ -2806,6 +2810,12 @@ 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 @@ -2912,6 +2922,7 @@ 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) @@ -3107,7 +3118,8 @@ endif enddo endif - endif + endif + call chainbuild_cart call kinetic(EK) write(iout,*) "just after kinetic" @@ -3371,8 +3383,11 @@ 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) @@ -3391,11 +3406,25 @@ 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 @@ -3427,20 +3456,26 @@ 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 @@ -3459,7 +3494,11 @@ 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 @@ -4436,7 +4475,7 @@ 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)) @@ -4452,6 +4491,7 @@ 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) @@ -4470,6 +4510,7 @@ 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) @@ -4488,9 +4529,12 @@ 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) @@ -4504,7 +4548,9 @@ 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) @@ -4751,12 +4797,14 @@ 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) @@ -4772,7 +4820,7 @@ 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) @@ -4943,8 +4991,8 @@ 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 @@ -4976,7 +5024,7 @@ 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 @@ -5029,8 +5077,8 @@ 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 @@ -5930,10 +5978,11 @@ ! 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 @@ -5948,6 +5997,17 @@ 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 @@ -6175,6 +6235,7 @@ endif #endif return +#undef DEBUG end subroutine friction_force !----------------------------------------------------------------------------- subroutine setup_fricmat @@ -6324,7 +6385,8 @@ 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 @@ -6601,10 +6663,11 @@ #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 @@ -6616,6 +6679,9 @@ #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 @@ -6631,6 +6697,14 @@ 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 @@ -6651,7 +6725,7 @@ 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 @@ -6811,6 +6885,7 @@ endif #endif return +#undef DEBUG end subroutine stochastic_force !----------------------------------------------------------------------------- subroutine sdarea(gamvec) @@ -7477,7 +7552,7 @@ 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 diff --git a/source/unres/REMD.F90 b/source/unres/REMD.F90 index 0405dbb..3cfe7dd 100644 --- a/source/unres/REMD.F90 +++ b/source/unres/REMD.F90 @@ -391,7 +391,7 @@ 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) @@ -1095,13 +1095,14 @@ 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 @@ -1320,10 +1321,10 @@ 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) diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index 1dfe01e..941735d 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -173,7 +173,7 @@ 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)) @@ -330,7 +330,7 @@ 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) @@ -905,7 +905,8 @@ 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 diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 599a108..80b281b 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -549,6 +549,11 @@ 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 @@ -1154,17 +1159,17 @@ -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 @@ -1961,6 +1966,7 @@ ! 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,& @@ -1968,6 +1974,9 @@ #else link_start=1 link_end=nhpb + clink_start=1 + clink_end=cnhpb + #endif return end subroutine hpb_partition diff --git a/source/unres/data/MD_data.F90 b/source/unres/data/MD_data.F90 index 898d55e..7959463 100644 --- a/source/unres/data/MD_data.F90 +++ b/source/unres/data/MD_data.F90 @@ -106,7 +106,7 @@ !----------------------------------------------------------------------------- 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 diff --git a/source/unres/data/control_data.F90 b/source/unres/data/control_data.F90 index b1e824d..e705b14 100644 --- a/source/unres/data/control_data.F90 +++ b/source/unres/data/control_data.F90 @@ -30,7 +30,8 @@ ! 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,& diff --git a/source/unres/data/energy_data.F90 b/source/unres/data/energy_data.F90 index fb4ff11..f09ac44 100644 --- a/source/unres/data/energy_data.F90 +++ b/source/unres/data/energy_data.F90 @@ -249,12 +249,14 @@ 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) diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index bb7d08d..c14e4ad 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -418,8 +418,9 @@ 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 @@ -430,6 +431,7 @@ 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 @@ -440,6 +442,7 @@ ! 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" @@ -5759,13 +5762,13 @@ 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 @@ -5801,22 +5804,22 @@ 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 @@ -5841,16 +5844,16 @@ 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 @@ -5869,13 +5872,14 @@ 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 @@ -12015,10 +12019,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! #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 @@ -12029,7 +12033,15 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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' @@ -12276,6 +12288,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! enddo ! enddo !el#define DEBUG +! print *,"gradbufc after summing" #ifdef DEBUG write (iout,*) "gradbufc after summing" do i=1,nres @@ -26950,6 +26963,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! !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 @@ -29764,7 +29782,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -29869,13 +29887,14 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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)) @@ -29997,7 +30016,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -30111,15 +30130,13 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -30136,6 +30153,17 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -30567,8 +30595,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -30579,7 +30607,10 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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 @@ -32573,12 +32604,19 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -32591,19 +32629,26 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -32919,7 +32964,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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 @@ -33046,11 +33091,11 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -33132,6 +33177,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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) @@ -33151,6 +33198,8 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 @@ -33186,6 +33235,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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)') & @@ -33204,6 +33254,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! ! 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 @@ -34253,7 +34304,7 @@ C !!!!!!!! NIE CZYTANE !!!!!!!!!!! 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 diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 index 866f5e2..4d74835 100644 --- a/source/unres/geometry.F90 +++ b/source/unres/geometry.F90 @@ -3828,8 +3828,48 @@ 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 @@ -3887,10 +3927,10 @@ 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)) @@ -3907,7 +3947,7 @@ 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 @@ -3974,6 +4014,21 @@ ! 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 diff --git a/source/unres/io.F90 b/source/unres/io.F90 index 3cb7f1c..64b4ede 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -1394,7 +1394,12 @@ 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 @@ -2678,8 +2683,8 @@ 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 @@ -2735,9 +2740,11 @@ 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 diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90 index f6058d7..a77c487 100644 --- a/source/unres/io_base.F90 +++ b/source/unres/io_base.F90 @@ -530,17 +530,18 @@ ! 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) @@ -587,33 +588,34 @@ 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 @@ -630,18 +632,18 @@ 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 @@ -653,33 +655,33 @@ ! 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 @@ -701,21 +703,21 @@ 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 @@ -723,11 +725,11 @@ 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 diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index cdd6231..42c794f 100644 --- a/source/unres/io_config.F90 +++ b/source/unres/io_config.F90 @@ -2761,8 +2761,14 @@ 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 @@ -4547,7 +4553,13 @@ ! 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" @@ -4559,7 +4571,7 @@ (c(j,ires+nres),j=1,3) enddo endif - + print *,seqalingbegin,nres if(.not.allocated(vbld)) then allocate(vbld(2*nres)) @@ -4864,6 +4876,7 @@ 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) diff --git a/source/unres/math.F90 b/source/unres/math.F90 index 798f90f..f6a800a 100644 --- a/source/unres/math.F90 +++ b/source/unres/math.F90 @@ -28,6 +28,7 @@ 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 @@ -101,6 +102,7 @@ 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 @@ -112,6 +114,7 @@ 17 AIIMIN=AJJ(I) IT=I 16 CONTINUE + print *,IT,"IT" IN=L AII(IN)=AIIMIN AJJ(IT)=1.0E+30 diff --git a/source/unres/md_calc.F90 b/source/unres/md_calc.F90 index 8b8e04d..ca9c526 100644 --- a/source/unres/md_calc.F90 +++ b/source/unres/md_calc.F90 @@ -1647,7 +1647,7 @@ 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 @@ -1676,6 +1676,7 @@ ! 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 ----- diff --git a/source/unres/minim.F90 b/source/unres/minim.F90 index f15d802..c661960 100644 --- a/source/unres/minim.F90 +++ b/source/unres/minim.F90 @@ -3662,13 +3662,13 @@ 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 @@ -3724,16 +3724,27 @@ 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 @@ -3742,7 +3753,8 @@ 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) @@ -3751,6 +3763,7 @@ enddo nvarx=k call chainbuild_cart +! print *,"before etotal" call etotal(energia(0)) call enerprint(energia(0)) #ifdef LBFGS @@ -3776,6 +3789,7 @@ 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 @@ -3838,10 +3852,11 @@ 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 @@ -3861,7 +3876,7 @@ 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 @@ -6800,6 +6815,7 @@ !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 @@ -6833,6 +6849,7 @@ 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 @@ -6919,8 +6936,10 @@ 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 diff --git a/source/unres/search.F90 b/source/unres/search.F90 index d85a192..4439b32 100644 --- a/source/unres/search.F90 +++ b/source/unres/search.F90 @@ -115,6 +115,7 @@ ! ! 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 diff --git a/source/xdrfpdb/xdrf2pdb-m.F90 b/source/xdrfpdb/xdrf2pdb-m.F90 index fe6eea7..b31e814 100644 --- a/source/xdrfpdb/xdrf2pdb-m.F90 +++ b/source/xdrfpdb/xdrf2pdb-m.F90 @@ -32,6 +32,9 @@ 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",&