X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=e34b3840191dabba240fbe66f40af5d196af8eb6;hb=10689ab7d813dfbdbb0c6e631d90234d78ea306a;hp=a436a938fcb5e92345acc3e4d9576f4730d33bcd;hpb=35f220f409bd5d21be33a402d79da2c23d3e0c3a;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index a436a93..e34b384 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -24,8 +24,18 @@ d_t_work_new,d_af_work,d_as_work,kinetic_force !(MAXRES6) real(kind=8),dimension(:,:),allocatable :: d_t_new,& d_a_old,d_a_short!,d_a !(3,0:MAXRES2) +! real(kind=8),dimension(:),allocatable :: d_a_work !(6*MAXRES) +! real(kind=8),dimension(:,:),allocatable :: Gmat,Ginv,A,& +! Gsqrp,Gsqrm,Gvec !(maxres2,maxres2) +! real(kind=8),dimension(:),allocatable :: Geigen !(maxres2) +! integer :: dimen,dimen1,dimen3 +! integer :: lang,count_reset_moment,count_reset_vel +! logical :: reset_moment,reset_vel,rattle,RESPA ! common /inertia/ ! common /langevin/ +! real(kind=8) :: rwat,etawat,stdfp,cPoise +! real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1) +! real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp) real(kind=8),dimension(:),allocatable :: stdforcp,stdforcsc !(MAXRES) !----------------------------------------------------------------------------- ! 'sizes.i' @@ -88,33 +98,33 @@ !el integer maxamino,maxnuc,maxbnd !el integer maxang,maxtors,maxpi !el integer maxpib,maxpit - integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres - integer,parameter :: maxval=8 - integer,parameter :: maxgrp=1000 - integer,parameter :: maxtyp=3000 - integer,parameter :: maxclass=500 - integer,parameter :: maxkey=10000 - integer,parameter :: maxrot=1000 - integer,parameter :: maxopt=1000 - integer,parameter :: maxhess=1000000 - integer :: maxlight !=8*maxatm - integer,parameter :: maxvib=1000 - integer,parameter :: maxgeo=1000 - integer,parameter :: maxcell=10000 - integer,parameter :: maxring=10000 - integer,parameter :: maxfix=10000 - integer,parameter :: maxbio=10000 - integer,parameter :: maxamino=31 - integer,parameter :: maxnuc=12 - integer :: maxbnd !=2*maxatm - integer :: maxang !=3*maxatm - integer :: maxtors !=4*maxatm - integer,parameter :: maxpi=100 - integer,parameter :: maxpib=2*maxpi - integer,parameter :: maxpit=4*maxpi +! integer :: maxatm !=2*nres !maxres2 maxres2=2*maxres +! integer,parameter :: maxval=8 +! integer,parameter :: maxgrp=1000 +! integer,parameter :: maxtyp=3000 +! integer,parameter :: maxclass=500 +! integer,parameter :: maxkey=10000 +! integer,parameter :: maxrot=1000 +! integer,parameter :: maxopt=1000 +! integer,parameter :: maxhess=1000000 +! integer :: maxlight !=8*maxatm +! integer,parameter :: maxvib=1000 +! integer,parameter :: maxgeo=1000 +! integer,parameter :: maxcell=10000 +! integer,parameter :: maxring=10000 +! integer,parameter :: maxfix=10000 +! integer,parameter :: maxbio=10000 +! integer,parameter :: maxamino=31 +! integer,parameter :: maxnuc=12 +! integer :: maxbnd !=2*maxatm +! integer :: maxang !=3*maxatm +! integer :: maxtors !=4*maxatm +! integer,parameter :: maxpi=100 +! integer,parameter :: maxpib=2*maxpi +! integer,parameter :: maxpit=4*maxpi !----------------------------------------------------------------------------- ! Maximum number of seed - integer,parameter :: max_seed=1 +! integer,parameter :: max_seed=1 !----------------------------------------------------------------------------- real(kind=8),dimension(:),allocatable :: stochforcvec !(MAXRES6) maxres6=6*maxres ! common /stochcalc/ stochforcvec @@ -131,7 +141,7 @@ real(kind=8),dimension(:,:),allocatable :: ginvfric !(2*nres,2*nres) !maxres2=2*maxres !----------------------------------------------------------------------------- ! common /przechowalnia/ subroutine: setup_fricmat -!el real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) + real(kind=8),dimension(:,:),allocatable :: fcopy !(2*nres,2*nres) !----------------------------------------------------------------------------- ! ! @@ -149,6 +159,7 @@ use control, only: tcpu use control_data use energy_data +! use io_conf, only:cartprint ! include 'DIMENSIONS' #ifdef MPI include 'mpif.h' @@ -176,6 +187,9 @@ real(kind=8),dimension(6*nres,2*nres) :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) real(kind=8),dimension(2*nres,2*nres) :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres real(kind=8),dimension(6*nres,6*nres) :: Pmat !(maxres6,maxres6) maxres6=6*maxres +! real(kind=8),dimension(:,:),allocatable :: Bmat,GBmat,Tmat !(MAXRES6,MAXRES2) (maxres2=2*maxres,maxres6=6*maxres) +! real(kind=8),dimension(:,:),allocatable :: Cmat_,Cinv !(maxres2,maxres2) maxres2=2*maxres +! real(kind=8),dimension(:,:),allocatable :: Pmat !(maxres6,maxres6) maxres6=6*maxres real(kind=8),dimension(6*nres) :: Td !(maxres6) maxres6=6*maxres real(kind=8),dimension(2*nres) :: ppvec !(maxres2) maxres2=2*maxres !el common /stochcalc/ stochforcvec @@ -191,12 +205,18 @@ logical :: osob nres2=2*nres nres6=6*nres +! if (.not.allocated(Bmat)) allocate(Bmat(nres6,nres2)) +! if (.not.allocated(GBmat)) allocate (GBmat(nres6,nres2)) +! if (.not.allocated(Tmat)) allocate (Tmat(nres6,nres2)) +! if (.not.allocated(Cmat_)) allocate(Cmat_(nres2,nres2)) +! if (.not.allocated(Cinv)) allocate (Cinv(nres2,nres2)) +! if (.not.allocated(Pmat)) allocate(Pmat(6*nres,6*nres)) if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres nbond=nct-nnt do i=nnt,nct - if (itype(i).ne.10) nbond=nbond+1 + if (itype(i,1).ne.10) nbond=nbond+1 enddo ! if (lprn1) then @@ -218,7 +238,7 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 Bmat(ind+j,ind1)=dC_norm(j,i+nres) @@ -295,9 +315,9 @@ Td(i)=Td(i)+vbl*Tmat(i,ind) enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then ind=ind+1 - Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind) + Td(i)=Td(i)+vbldsc0(1,itype(k,1))*Tmat(i,ind) endif enddo enddo @@ -328,7 +348,7 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then do j=1,3 ind=ind+1 zapas(ind)=-gxcart(j,i)+stochforcvec(ind) @@ -399,9 +419,9 @@ i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 - xx=vbld(i+nres)-vbldsc0(1,itype(i)) + xx=vbld(i+nres)-vbldsc0(1,itype(i,1)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif @@ -427,11 +447,11 @@ endif enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 blen2 = scalar(dc(1,i+nres),dc(1,i+nres)) - ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2 - diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2)) + ppvec(ind)=2*vbldsc0(1,itype(i,1))**2-blen2 + diffbond=dabs(vbldsc0(1,itype(i,1))-dsqrt(blen2)) if (diffbond.gt.diffmax) diffmax=diffbond if (ppvec(ind).gt.0.0d0) then ppvec(ind)=dsqrt(ppvec(ind)) @@ -470,7 +490,7 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then do j=1,3 dc(j,i+nres)=zapas(ind+j) dc_work(ind+j)=zapas(ind+j) @@ -514,9 +534,9 @@ i,(dC(j,i),j=1,3),xx enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 - xx=vbld(i+nres)-vbldsc0(1,itype(i)) + xx=vbld(i+nres)-vbldsc0(1,itype(i,1)) write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),xx endif @@ -537,6 +557,7 @@ potE=potEcomp(0)-potEcomp(20) call cartgrad totT=totT+d_time + totTafm=totT ! Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) #ifdef MPI @@ -640,26 +661,34 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.MD' ! include 'COMMON.IOUNITS' - real(kind=8) :: KE_total + real(kind=8) :: KE_total,mscab - integer :: i,j,k,iti + integer :: i,j,k,iti,mnum real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),& mag1,mag2,v(3) - +#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 KEt_p=0.0d0 KEt_sc=0.0d0 -! write (iout,*) "ISC",(isc(itype(i)),i=1,nres) +! write (iout,*) "ISC",(isc(itype(i,1)),i=1,nres) ! The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) ! write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) - KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + KEt_p=KEt_p+mp(mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo @@ -672,8 +701,16 @@ incr(j)=d_t(j,0) enddo do i=nnt,nct - iti=iabs(itype(i)) - if (itype(i).eq.10) then + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (mnum.eq.5) then + mscab=0.0d0 + else + mscab=msc(iti,mnum) + endif +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.(mnum.eq.5)) then do j=1,3 v(j)=incr(j) enddo @@ -684,7 +721,7 @@ endif ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) - KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(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) do j=1,3 incr(j)=incr(j)+d_t(j,i) @@ -695,13 +732,14 @@ ! The part due to stretching and rotation of the peptide groups KEr_p=0.0D0 do i=nnt,nct-1 + mnum=molnum(i) ! write (iout,*) "i",i ! write (iout,*) "i",i," mag1",mag1," mag2",mag2 do j=1,3 incr(j)=d_t(j,i) enddo ! write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) - KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) & + KEr_p=KEr_p+Ip(mnum)*(incr(1)*incr(1)+incr(2)*incr(2) & +incr(3)*incr(3)) enddo ! goto 111 @@ -709,20 +747,23 @@ ! The rotational part of the side chain virtual bond KEr_sc=0.0D0 do i=nnt,nct - iti=iabs(itype(i)) - if (itype(i).ne.10) then + mnum=molnum(i) + iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 incr(j)=d_t(j,nres+i) enddo ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) - KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ & + 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 - KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc) + KE_total=0.5d0*(KEt_p+KEt_sc+0.25d0*KEr_p+KEr_sc) ! write (iout,*) "KE_total",KE_total return end subroutine kinetic @@ -734,7 +775,9 @@ ! The driver for molecular dynamics subroutines !------------------------------------------------ use comm_gucio +! use MPI use control, only:tcpu,ovrtim +! use io_comm, only:ilen use control_data use compare, only:secondary2,hairpin use io, only:cartout,statout @@ -776,8 +819,10 @@ real(kind=8) :: tt0,scalfac integer :: nres2 nres2=2*nres + print *, "ENTER MD" ! #ifdef MPI + print *,"MY tmpdir",tmpdir,ilen(tmpdir) if (ilen(tmpdir).gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//"_" & //liczba(:ilen(liczba))//'.rst') @@ -796,10 +841,14 @@ #else tt0 = tcpu() #endif + print *,"just befor setup matix",nres ! Determine the inverse of the inertia matrix. call setup_MD_matrices ! Initialize MD + print *,"AFTER SETUP MATRICES" call init_MD + print *,"AFTER INIT MD" + #ifdef MPI t_MDsetup = MPI_Wtime()-tt0 #else @@ -843,7 +892,9 @@ stop #endif else if (lang.eq.1 .or. lang.eq.4) then + print *,"before setup_fricmat" call setup_fricmat + print *,"after setup_fricmat" endif #ifdef MPI t_langsetup=MPI_Wtime()-tt0 @@ -936,7 +987,11 @@ #endif endif if (ntwe.ne.0) then - if (mod(itime,ntwe).eq.0) call statout(itime) + if (mod(itime,ntwe).eq.0) then + call statout(itime) + call returnbox +! call check_ecartint + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -949,7 +1004,8 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) @@ -970,6 +1026,7 @@ #endif endif if (mod(itime,ntwx).eq.0) then + call returnbox write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) @@ -982,10 +1039,13 @@ if (rstcount.eq.1000.or.itime.eq.n_timestep) then open(irest2,file=rest2name,status='unknown') write(irest2,*) totT,EK,potE,totE,t_bath - do i=1,2*nres + totTafm=totT +! AL 4/17/17: Now writing d_t(0,:) too + do i=0,2*nres write (irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres +! AL 4/17/17: Now writing d_c(0,:) too + do i=0,2*nres write (irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo close(irest2) @@ -1142,6 +1202,12 @@ ! Calculate energy and forces call zerograd call etotal(potEcomp) +! AL 4/17/17: Reduce the steps if NaNs occurred. + if (potEcomp(0).gt.0.99e20 .or. isnan(potEcomp(0)).gt.0) then + d_time=d_time/2 + cycle + endif +! end change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(potEcomp) #ifdef TIMING_ENE @@ -1260,6 +1326,7 @@ endif if (rattle) call rattle2 totT=totT+d_time + totTafm=totT if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -1324,8 +1391,10 @@ ! include 'DIMENSIONS' use comm_gucio use comm_cipiszcze +! use MPI use control, only:tcpu use control_data +! use io_conf, only:cartprint #ifdef MPI include 'mpif.h' integer :: IERROR,ERRCODE @@ -1515,6 +1584,25 @@ ! Calculate energy and forces call zerograd call etotal_short(energia_short) +! AL 4/17/17: Exit itime_split loop when energy goes infinite + if (energia_short(0).gt.0.99e20 .or. isnan(energia_short(0)) ) then + if (PRINT_AMTS_MSG) & + write (iout,*) "Infinities/NaNs in energia_short",energia_short(0),"; increasing ntime_split to",ntime_split + ntime_split=ntime_split*2 + if (ntime_split.gt.maxtime_split) then +#ifdef MPI + write (iout,*) & + "Cannot rescue the run; aborting job. Retry with a smaller time step" + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#else + write (iout,*) & + "Cannot rescue the run; terminating. Retry with a smaller time step" +#endif + endif + exit + endif +! End change if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_short) #ifdef TIMING_ENE @@ -1557,6 +1645,8 @@ if (ntime_split.lt.maxtime_split) then scale=.true. ntime_split=ntime_split*2 +! AL 4/17/17: We should exit the itime_split loop when acceleration change is too big + exit do i=0,2*nres do j=1,3 dc_old(j,i)=dc_old0(j,i) @@ -1619,6 +1709,17 @@ #endif call zerograd call etotal_long(energia_long) + if (energia_long(0).gt.0.99e20 .or. isnan(energia_long(0))) then +#ifdef MPI + write (iout,*) & + "Infinitied/NaNs in energia_long, Aborting MPI job." + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE) +#else + write (iout,*) "Infinitied/NaNs in energia_long, terminating." + stop +#endif + endif if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_long) #ifdef TIMING_ENE @@ -1665,6 +1766,7 @@ potE=potEcomp(0)-potEcomp(20) ! potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT ! Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -1702,7 +1804,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - integer :: i,j,inres + integer :: i,j,inres,mnum do j=1,3 d_t(j,0)=d_t(j,0)+0.5d0*d_a(j,0)*d_time @@ -1713,7 +1815,10 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1739,7 +1844,7 @@ ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8) :: adt,adt2 - integer :: i,j,inres + integer :: i,j,inres,mnum #ifdef DEBUG write (iout,*) "VELVERLET1 START: DC" @@ -1765,7 +1870,10 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time @@ -1801,7 +1909,7 @@ ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - integer :: i,j,inres + integer :: i,j,inres,mnum do j=1,3 d_t(j,0)=d_t_new(j,0)+0.5d0*d_a(j,0)*d_time @@ -1812,7 +1920,11 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1895,7 +2007,7 @@ ! position and velocity increments included. real(kind=8) :: sqrt13 = 0.57735026918962576451d0 ! 1/sqrt(3) real(kind=8) :: adt,adt2 - integer :: i,j,ind,inres + integer :: i,j,ind,inres,mnum ! ! Add the contribution from BOTH friction and stochastic force to the ! coordinates, but ONLY the contribution from the friction forces to velocities @@ -1919,7 +2031,11 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time @@ -1957,7 +2073,7 @@ ! include 'COMMON.NAMES' real(kind=8),dimension(6*nres) :: stochforcvec,d_as_work1 !(MAXRES6) maxres6=6*maxres real(kind=8) :: cos60 = 0.5d0, sin60 = 0.86602540378443864676d0 - integer :: i,j,ind,inres + integer :: i,j,ind,inres,mnum ! Revised 3/31/05 AL: correlation between random contributions to ! position and velocity increments included. ! The correlation coefficients are calculated at low-friction limit. @@ -1987,7 +2103,11 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) & @@ -2019,7 +2139,7 @@ ! include 'COMMON.IOUNITS' real(kind=8),dimension(3) :: aux,accel,accel_old real(kind=8) :: dacc - integer :: i,j + integer :: i,j,mnum do j=1,3 ! aux(j)=d_a(j,0)-d_a_old(j,0) @@ -2057,7 +2177,11 @@ enddo endif do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 ! accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) accel_old(j)=accel_old(j)+d_a_old(j,i+nres) @@ -2114,7 +2238,8 @@ enddo endif ! Side chains - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.& + molnum(i).ne.5) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) @@ -2147,7 +2272,7 @@ ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' real(kind=8) :: T_half,fact - integer :: i,j,inres + integer :: i,j,inres,mnum ! T_half=2.0d0/(dimen3*Rb)*EK fact=dsqrt(1.0d0+(d_time/tau_bath)*(t_bath/T_half-1.0d0)) @@ -2163,7 +2288,11 @@ enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) @@ -2218,7 +2347,7 @@ character(len=50) :: tytul logical :: file_exist !el common /gucio/ cm - integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr + integer :: i,j,ipos,iq,iw,nft_sc,iretcode,nfun,itime,ierr,mnum real(kind=8) :: etot,tt0 logical :: fail @@ -2228,11 +2357,13 @@ ! if the friction coefficients do not depend on surface area if (lang.gt.0 .and. .not.surfarea) then do i=nnt,nct-1 - stdforcp(i)=stdfp*dsqrt(gamp) + mnum=molnum(i) + stdforcp(i)=stdfp(mnum)*dsqrt(gamp(mnum)) enddo do i=nnt,nct - stdforcsc(i)=stdfsc(iabs(itype(i))) & - *dsqrt(gamsc(iabs(itype(i)))) + mnum=molnum(i) + stdforcsc(i)=stdfsc(iabs(itype(i,mnum)),mnum) & + *dsqrt(gamsc(iabs(itype(i,mnum)),mnum)) enddo endif ! Open the pdb file for snapshotshots @@ -2342,6 +2473,7 @@ endif call random_vel totT=0.0d0 + totTafm=totT endif else ! Generate initial velocities @@ -2349,6 +2481,7 @@ write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 + totTafm=totT endif ! rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then @@ -2375,8 +2508,8 @@ write (iout,*) (vcm(j),j=1,3) endif if (.not.rest) then - call chainbuild - if(iranconf.ne.0) then +! call chainbuild + if(iranconf.ne.0 .or.indpdb.gt.0.and..not.unres_pdb .or.preminim) then if (overlapsc) then print *, 'Calling OVERLAP_SC' call overlap_sc(fail) @@ -2402,14 +2535,11 @@ endif endif call chainbuild_cart -write(iout,*) "przed kinetic, EK=",EK call kinetic(EK) if (tbf) then call verlet_bath endif -write(iout,*) "dimen3",dimen3,"Rb",Rb,"EK",EK kinetic_T=2.0d0/(dimen3*Rb)*EK -write(iout,*) "kinetic_T",kinetic_T if(me.eq.king.or..not.out1file)then call cartprint call intout @@ -2431,9 +2561,7 @@ write(iout,*) "kinetic_T",kinetic_T #endif potE=potEcomp(0) call cartgrad -write(iout,*) "kinetic_T if large",kinetic_T call lagrangian -write(iout,*) "kinetic_T if large",kinetic_T call max_accel if (amax*d_time .gt. dvmax) then d_time=d_time*dvmax/amax @@ -2504,9 +2632,7 @@ write(iout,*) "kinetic_T if large",kinetic_T #endif #endif call cartgrad -write(iout,*) "przed lagrangian" call lagrangian -write(iout,*) "po lagrangian" if(.not.out1file .and. large) then write (iout,*) "energia_long",energia_long(0),& " energia_short",energia_short(0),& @@ -2539,9 +2665,7 @@ write(iout,*) "po lagrangian" #endif #endif call cartgrad -write(iout,*) "przed lagrangian2" call lagrangian -write(iout,*) "po lagrangian2" if(.not.out1file .and. large) then write (iout,*) "energia_long",energia_long(0) write (iout,*) "Initial slow-force accelerations" @@ -2556,7 +2680,6 @@ write(iout,*) "po lagrangian2" t_enegrad=t_enegrad+tcpu()-tt0 #endif endif -write(iout,*) "end init MD" return end subroutine init_MD !----------------------------------------------------------------------------- @@ -2565,6 +2688,7 @@ write(iout,*) "end init MD" ! implicit real*8 (a-h,o-z) use energy_data use random, only:anorm_distr + use MD_data ! include 'DIMENSIONS' ! include 'COMMON.CONTROL' ! include 'COMMON.VAR' @@ -2583,7 +2707,17 @@ write(iout,*) "end init MD" ! include 'COMMON.NAMES' ! include 'COMMON.TIME1' real(kind=8) :: xv,sigv,lowb,highb ,Ek1 - integer :: i,j,ii,k,ind +!#define DEBUG +#ifdef FIVEDIAG + real(kind=8) ,allocatable, dimension(:) :: DDU1,DDU2,DL2,DL1,xsolv,DML,rs + real(kind=8) :: sumx +#ifdef DEBUG + real(kind=8) ,allocatable, dimension(:) :: rsold + real (kind=8),allocatable,dimension(:,:) :: matold + integer :: iti +#endif +#endif + integer :: i,j,ii,k,ind,mark,imark,mnum ! Generate random velocities from Gaussian distribution of mean 0 and std of KT/m ! First generate velocities in the eigenspace of the G matrix ! write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 @@ -2597,22 +2731,317 @@ write(iout,*) "end init MD" lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) -! write (iout,*) "i",i," ii",ii," geigen",geigen(i),& -! " d_t_work_new",d_t_work_new(ii) +#ifdef DEBUG + write (iout,*) "i",i," ii",ii," geigen",geigen(i),& + " d_t_work_new",d_t_work_new(ii) +#endif enddo enddo +#ifdef DEBUG ! diagnostics -! Ek1=0.0d0 -! ii=0 -! do i=1,dimen -! do k=1,3 -! ii=ii+1 -! Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 -! enddo -! enddo -! write (iout,*) "Ek from eigenvectors",Ek1 + Ek1=0.0d0 + ii=0 + do i=1,dimen + do k=1,3 + ii=ii+1 + Ek1=Ek1+0.5d0*geigen(i)*d_t_work_new(ii)**2 + enddo + enddo + write (iout,*) "Ek from eigenvectors",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) ! end diagnostics +#endif +#ifdef FIVEDIAG ! Transform velocities to UNRES coordinate space + allocate (DL1(2*nres)) + allocate (DDU1(2*nres)) + allocate (DL2(2*nres)) + allocate (DDU2(2*nres)) + allocate (xsolv(2*nres)) + allocate (DML(2*nres)) + allocate (rs(2*nres)) +#ifdef DEBUG + allocate (rsold(2*nres)) + allocate (matold(2*nres,2*nres)) + matold=0 + matold(1,1)=DMorig(1) + matold(1,2)=DU1orig(1) + matold(1,3)=DU2orig(1) + write (*,*) DMorig(1),DU1orig(1),DU2orig(1) + + do i=2,dimen-1 + do j=1,dimen + if (i.eq.j) then + matold(i,j)=DMorig(i) + matold(i,j-1)=DU1orig(i-1) + if (j.ge.3) then + matold(i,j-2)=DU2orig(i-2) + endif + + endif + + enddo + do j=1,dimen-1 + if (i.eq.j) then + matold(i,j+1)=DU1orig(i) + + end if + enddo + do j=1,dimen-2 + if(i.eq.j) then + matold(i,j+2)=DU2orig(i) + endif + enddo + enddo + matold(dimen,dimen)=DMorig(dimen) + matold(dimen,dimen-1)=DU1orig(dimen-1) + matold(dimen,dimen-2)=DU2orig(dimen-2) + write(iout,*) "old gmatrix" + call matout(dimen,dimen,2*nres,2*nres,matold) +#endif + d_t_work=0.0d0 + do i=1,dimen +! Find the ith eigenvector of the pentadiagonal inertiq matrix + do imark=dimen,1,-1 + + do j=1,imark-1 + DML(j)=DMorig(j)-geigen(i) + enddo + do j=imark+1,dimen + DML(j-1)=DMorig(j)-geigen(i) + enddo + do j=1,imark-2 + DDU1(j)=DU1orig(j) + enddo + DDU1(imark-1)=DU2orig(imark-1) + do j=imark+1,dimen-1 + DDU1(j-1)=DU1orig(j) + enddo + do j=1,imark-3 + DDU2(j)=DU2orig(j) + enddo + DDU2(imark-2)=0.0d0 + DDU2(imark-1)=0.0d0 + do j=imark,dimen-3 + DDU2(j)=DU2orig(j+1) + enddo + do j=1,dimen-3 + DL2(j+2)=DDU2(j) + enddo + do j=1,dimen-2 + DL1(j+1)=DDU1(j) + enddo +#ifdef DEBUG + write (iout,*) "DL2,DL1,DML,DDU1,DDU2" + write(iout,'(10f10.5)') (DL2(k),k=3,dimen-1) + write(iout,'(10f10.5)') (DL1(k),k=2,dimen-1) + write(iout,'(10f10.5)') (DML(k),k=1,dimen-1) + write(iout,'(10f10.5)') (DDU1(k),k=1,dimen-2) + write(iout,'(10f10.5)') (DDU2(k),k=1,dimen-3) +#endif + rs=0.0d0 + if (imark.gt.2) rs(imark-2)=-DU2orig(imark-2) + if (imark.gt.1) rs(imark-1)=-DU1orig(imark-1) + if (imark.lt.dimen) rs(imark)=-DU1orig(imark) + if (imark.lt.dimen-1) rs(imark+1)=-DU2orig(imark) +#ifdef DEBUG + rsold=rs +#endif +! write (iout,*) "Vector rs" +! do j=1,dimen-1 +! write (iout,*) j,rs(j) +! enddo + xsolv=0.0d0 + call FDIAG(dimen-1,DL2,DL1,DML,DDU1,DDU2,rs,xsolv,mark) + + if (mark.eq.1) then + +#ifdef DEBUG +! Check solution + do j=1,imark-1 + sumx=-geigen(i)*xsolv(j) + do k=1,imark-1 + sumx=sumx+matold(j,k)*xsolv(k) + enddo + do k=imark+1,dimen + sumx=sumx+matold(j,k)*xsolv(k-1) + enddo + write(iout,'(i5,3f10.5)') j,sumx,rsold(j),sumx-rsold(j) + enddo + do j=imark+1,dimen + sumx=-geigen(i)*xsolv(j-1) + do k=1,imark-1 + sumx=sumx+matold(j,k)*xsolv(k) + enddo + do k=imark+1,dimen + sumx=sumx+matold(j,k)*xsolv(k-1) + enddo + write(iout,'(i5,3f10.5)') j-1,sumx,rsold(j-1),sumx-rsold(j-1) + enddo +! end check + write (iout,*)& + "Solution of equations system",i," for eigenvalue",geigen(i) + do j=1,dimen-1 + write(iout,'(i5,f10.5)') j,xsolv(j) + enddo +#endif + do j=dimen-1,imark,-1 + xsolv(j+1)=xsolv(j) + enddo + xsolv(imark)=1.0d0 +#ifdef DEBUG + write (iout,*) "Un-normalized eigenvector",i," for eigenvalue",geigen(i) + do j=1,dimen + write(iout,'(i5,f10.5)') j,xsolv(j) + enddo +#endif +! Normalize ith eigenvector + sumx=0.0d0 + do j=1,dimen + sumx=sumx+xsolv(j)**2 + enddo + sumx=dsqrt(sumx) + do j=1,dimen + xsolv(j)=xsolv(j)/sumx + enddo +#ifdef DEBUG + write (iout,*) "Eigenvector",i," for eigenvalue",geigen(i) + do j=1,dimen + write(iout,'(i5,3f10.5)') j,xsolv(j) + enddo +#endif +! All done at this point for eigenvector i, exit loop + exit + + endif ! mark.eq.1 + + enddo ! imark + + if (mark.ne.1) then + write (iout,*) "Unable to find eigenvector",i + endif + +! write (iout,*) "i=",i + do k=1,3 +! write (iout,*) "k=",k + do j=1,dimen + ind=(j-1)*3+k +! write(iout,*) "ind",ind," ind1",3*(i-1)+k + d_t_work(ind)=d_t_work(ind) & + +xsolv(j)*d_t_work_new(3*(i-1)+k) + enddo + enddo + enddo ! i (loop over the eigenvectors) + +#ifdef DEBUG + write (iout,*) "d_t_work" + do i=1,3*dimen + write (iout,"(i5,f10.5)") i,d_t_work(i) + enddo + Ek1=0.0d0 + ii=0 + do i=nnt,nct +! if (itype(i,1).eq.10) then + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + iti=iabs(itype(i,mnum)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.(mnum.eq.5)) then + j=ii+3 + else + j=ii+6 + endif + if (i.lt.nct) then + do k=1,3 +! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1 + Ek1=Ek1+0.5d0*mp(mnum)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& + 0.5d0*0.25d0*IP(mnum)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 + enddo + endif + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) ii=ii+3 + write (iout,*) "i",i," itype",itype(i,mnum)," mass",msc(itype(i,mnum),mnum) + write (iout,*) "ii",ii + do k=1,3 + ii=ii+1 + write (iout,*) "k",k," ii",ii,"EK1",EK1 + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5))& + Ek1=Ek1+0.5d0*Isc(iabs(itype(i,mnum)),mnum)*(d_t_work(ii)-d_t_work(ii-3))**2 + Ek1=Ek1+0.5d0*msc(iabs(itype(i,mnum)),mnum)*d_t_work(ii)**2 + enddo + write (iout,*) "i",i," ii",ii + enddo + write (iout,*) "Ek from d_t_work",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) +#endif + deallocate(DDU1) + deallocate(DDU2) + deallocate(DL2) + deallocate(DL1) + deallocate(xsolv) + deallocate(DML) + deallocate(rs) +#ifdef DEBUG + deallocate(matold) + deallocate(rsold) +#endif + ind=1 + do j=nnt,nct + do k=1,3 + d_t(k,j)=d_t_work(ind) + ind=ind+1 + enddo + mnum=molnum(j) + if (itype(j,1).ne.10 .and. itype(j,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then + do k=1,3 + d_t(k,j+nres)=d_t_work(ind) + ind=ind+1 + enddo + end if + enddo +#ifdef DEBUG + write (iout,*) "Random velocities in the Calpha,SC space" + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + enddo + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) + enddo +#endif + do j=1,3 + d_t(j,0)=d_t(j,nnt) + enddo + do i=nnt,nct +! if (itype(i,1).eq.10) then + mnum=molnum(i) + if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)& + .or.(mnum.eq.5)) then + do j=1,3 + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + else + do j=1,3 + d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i) + d_t(j,i)=d_t(j,i+1)-d_t(j,i) + enddo + end if + enddo +#ifdef DEBUG + write (iout,*)"Random velocities in the virtual-bond-vector space" + do i=nnt,nct-1 + write (iout,'(i3,3f10.5)') i,(d_t(j,i),j=1,3) + enddo + do i=nnt,nct + write (iout,'(i3,3f10.5)') i,(d_t(j,i+nres),j=1,3) + enddo + call kinetic(Ek1) + write (iout,*) "Ek from d_t_work",Ek1 + write (iout,*) "Kinetic temperatures",2*Ek1/(3*dimen*Rb) +#endif +#else do k=0,2 do i=1,dimen ind=(i-1)*3+k+1 @@ -2637,17 +3066,22 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) enddo endif enddo +#endif ! call kinetic(EK) ! write (iout,*) "Kinetic energy",Ek,EK1," kinetic temperature",& ! 2.0d0/(dimen3*Rb)*EK,2.0d0/(dimen3*Rb)*EK1 ! call flush(iout) +! write(iout,*) "end init MD" return end subroutine random_vel !----------------------------------------------------------------------------- @@ -2876,7 +3310,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -2924,7 +3361,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -2990,7 +3430,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -3154,7 +3597,7 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -3203,7 +3646,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -3269,7 +3715,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) +! if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -3305,13 +3754,13 @@ write(iout,*) "end init MD" real(kind=8),dimension(3,3) :: Im,Imcp,eigvec,Id real(kind=8),dimension(3) :: pr,eigval,L,vp,vrot - real(kind=8) :: M_SC,mag,mag2 + real(kind=8) :: M_SC,mag,mag2,M_PEP real(kind=8),dimension(3,0:nres) :: vpp !(3,0:MAXRES) real(kind=8),dimension(3) :: vs_p,pp,incr,v real(kind=8),dimension(3,3) :: pr1,pr2 !el common /gucio/ cm - integer :: iti,inres,i,j,k + integer :: iti,inres,i,j,k,mnum do i=1,3 do j=1,3 Im(i,j)=0.0d0 @@ -3323,91 +3772,106 @@ write(iout,*) "end init MD" vrot(i)=0.0d0 enddo ! calculating the center of the mass of the protein + M_PEP=0.0d0 do i=nnt,nct-1 + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle + M_PEP=M_PEP+mp(mnum) do j=1,3 - cm(j)=cm(j)+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*cm(j) - enddo +! do j=1,3 +! cm(j)=mp(1)*cm(j) +! enddo M_SC=0.0d0 do i=nnt,nct - iti=iabs(itype(i)) - M_SC=M_SC+msc(iabs(iti)) + mnum=molnum(i) + if (mnum.eq.5) cycle + iti=iabs(itype(i,mnum)) + M_SC=M_SC+msc(iabs(iti),mnum) inres=i+nres do j=1,3 - cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) + cm(j)=cm(j)+msc(iabs(iti),mnum)*c(j,inres) enddo enddo do j=1,3 - cm(j)=cm(j)/(M_SC+(nct-nnt)*mp) + cm(j)=cm(j)/(M_SC+M_PEP) enddo do i=nnt,nct-1 + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo - Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3)) - Im(1,2)=Im(1,2)-mp*pr(1)*pr(2) - Im(1,3)=Im(1,3)-mp*pr(1)*pr(3) - Im(2,3)=Im(2,3)-mp*pr(2)*pr(3) - Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1)) - Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2)) + Im(1,1)=Im(1,1)+mp(mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-mp(mnum)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-mp(mnum)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-mp(mnum)*pr(2)*pr(3) + 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 do i=nnt,nct - iti=iabs(itype(i)) + mnum=molnum(i) + iti=iabs(itype(i,mnum)) + if (mnum.eq.5) cycle inres=i+nres do j=1,3 pr(j)=c(j,inres)-cm(j) enddo - Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3)) - Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2) - Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3) - Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3) - Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1)) - Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2)) + Im(1,1)=Im(1,1)+msc(iabs(iti),mnum)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-msc(iabs(iti),mnum)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-msc(iabs(iti),mnum)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-msc(iabs(iti),mnum)*pr(2)*pr(3) + 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 do i=nnt,nct-1 - Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* & + mnum=molnum(i) + Im(1,1)=Im(1,1)+Ip(mnum)*(1-dc_norm(1,i)*dc_norm(1,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* & + Im(1,2)=Im(1,2)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* & + Im(1,3)=Im(1,3)+Ip(mnum)*(-dc_norm(1,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* & + Im(2,3)=Im(2,3)+Ip(mnum)*(-dc_norm(2,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* & + Im(2,2)=Im(2,2)+Ip(mnum)*(1-dc_norm(2,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* & + 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 do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then - iti=iabs(itype(i)) + mnum=molnum(i) +! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then + iti=iabs(itype(i,mnum)) inres=i+nres - Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* & + Im(1,1)=Im(1,1)+Isc(iti,mnum)*(1-dc_norm(1,inres)* & dc_norm(1,inres))*vbld(inres)*vbld(inres) - Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* & + Im(1,2)=Im(1,2)-Isc(iti,mnum)*(dc_norm(1,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) - Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* & + Im(1,3)=Im(1,3)-Isc(iti,mnum)*(dc_norm(1,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) - Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* & + Im(2,3)=Im(2,3)-Isc(iti,mnum)*(dc_norm(2,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) - Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* & + Im(2,2)=Im(2,2)+Isc(iti,mnum)*(1-dc_norm(2,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) - Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* & + Im(3,3)=Im(3,3)+Isc(iti,mnum)*(1-dc_norm(3,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo call angmom(cm,L) ! write(iout,*) "The angular momentum before adjustment:" -! write(iout,*) (L(j),j=1,3) +! write(iout,*) (L(j),j=1,3) Im(2,1)=Im(1,2) Im(3,1)=Im(1,3) @@ -3417,7 +3881,7 @@ write(iout,*) "end init MD" do i=1,3 do j=1,3 Imcp(i,j)=Im(i,j) - Id(i,j)=0.0d0 + Id(i,j)=0.0d0 enddo enddo @@ -3470,7 +3934,11 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then +! if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then +! if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then inres=i+nres call vecpr(vrot(1),dc(1,inres),vp) do j=1,3 @@ -3480,7 +3948,7 @@ write(iout,*) "end init MD" enddo call angmom(cm,L) ! write(iout,*) "The angular momentum after adjustment:" -! write(iout,*) (L(j),j=1,3) +! write(iout,*) (L(j),j=1,3) return end subroutine inertia_tensor @@ -3500,9 +3968,9 @@ write(iout,*) "end init MD" ! include 'COMMON.INTERACT' ! include 'COMMON.IOUNITS' ! include 'COMMON.NAMES' - + real(kind=8) :: mscab real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp - integer :: iti,inres,i,j + integer :: iti,inres,i,j,mnum ! Calculate the angular momentum do j=1,3 L(j)=0.0d0 @@ -3511,6 +3979,8 @@ write(iout,*) "end init MD" incr(j)=d_t(j,0) enddo do i=nnt,nct-1 + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo @@ -3522,7 +3992,7 @@ write(iout,*) "end init MD" enddo call vecpr(pr(1),v(1),vp) do j=1,3 - L(j)=L(j)+mp*vp(j) + L(j)=L(j)+mp(mnum)*vp(j) enddo do j=1,3 pr(j)=0.5d0*dc(j,i) @@ -3530,19 +4000,26 @@ write(iout,*) "end init MD" enddo call vecpr(pr(1),pp(1),vp) do j=1,3 - L(j)=L(j)+Ip*vp(j) + L(j)=L(j)+Ip(mnum)*vp(j) enddo enddo do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct - iti=iabs(itype(i)) + mnum=molnum(i) + iti=iabs(itype(i,mnum)) inres=i+nres + if (mnum.eq.5) then + mscab=0.0d0 + else + mscab=msc(iti,mnum) + endif do j=1,3 pr(j)=c(j,inres)-cm(j) enddo - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo @@ -3552,19 +4029,20 @@ write(iout,*) "end init MD" enddo endif 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)+msc(iabs(iti))*vp(j) + L(j)=L(j)+mscab*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 v(j)=incr(j)+d_t(j,inres) enddo call vecpr(dc(1,inres),d_t(1,inres),vp) do j=1,3 - L(j)=L(j)+Isc(iti)*vp(j) + L(j)=L(j)+Isc(iti,mnum)*vp(j) enddo endif do j=1,3 @@ -3589,7 +4067,7 @@ write(iout,*) "end init MD" ! include 'COMMON.IOUNITS' real(kind=8),dimension(3) :: vcm,vv real(kind=8) :: summas,amas - integer :: i,j + integer :: i,j,mnum do j=1,3 vcm(j)=0.0d0 @@ -3597,15 +4075,22 @@ write(iout,*) "end init MD" enddo summas=0.0d0 do i=nnt,nct + mnum=molnum(i) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) if (i.lt.nct) then - summas=summas+mp + summas=summas+mp(mnum) do j=1,3 - vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i)) + vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i)) enddo endif - amas=msc(iabs(itype(i))) + if (mnum.ne.5) then + amas=msc(iabs(itype(i,mnum)),mnum) + else + amas=0.0d0 + endif summas=summas+amas - if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then do j=1,3 vcm(j)=vcm(j)+amas*(vv(j)+d_t(j,i+nres)) enddo @@ -3674,7 +4159,9 @@ write(iout,*) "end init MD" if (lprn) write (iout,*) "RATTLE1" nbond=nct-nnt do i=nnt,nct - if (itype(i).ne.10) nbond=nbond+1 + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) nbond=nbond+1 enddo ! Make a folded form of the Ginv-matrix ind=0 @@ -3693,7 +4180,9 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + mnum=molnum(k) + if (itype(k,1).ne.10 .and. itype(k,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then jj=jj+1 do l=1,3 ind1=ind1+1 @@ -3704,7 +4193,9 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) ii=ii+1 do j=1,3 ind=ind+1 @@ -3718,7 +4209,7 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 @@ -3750,7 +4241,7 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i+nres) @@ -3766,7 +4257,7 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres) @@ -3781,9 +4272,9 @@ write(iout,*) "end init MD" x(ind)=vbld(i+1)**2-vbl**2 enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 - x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2 + x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 endif enddo if (lprn) then @@ -3800,7 +4291,10 @@ write(iout,*) "end init MD" i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then + ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& @@ -3854,7 +4348,7 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 do j=1,3 xx=0.0d0 @@ -3880,9 +4374,11 @@ write(iout,*) "end init MD" i,(dC(j,i),j=1,3),x(ind),xx enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) ind=ind+1 - xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2 + xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),x(ind),xx endif @@ -3895,7 +4391,7 @@ write(iout,*) "end init MD" i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& @@ -3970,7 +4466,9 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres) @@ -3999,7 +4497,9 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then ind=ind+1 do j=1,nbond Cmat(ind,j)=0.0d0 @@ -4016,7 +4516,9 @@ write(iout,*) "end init MD" x(ind)=scalar(d_t(1,i),dC(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then ind=ind+1 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres)) endif @@ -4030,7 +4532,9 @@ write(iout,*) "end init MD" i,ind,(d_t(j,i),j=1,3),x(ind) enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind) @@ -4065,7 +4569,9 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) then ind=ind+1 do j=1,3 xx=0.0d0 @@ -4086,7 +4592,9 @@ write(iout,*) "end init MD" i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)& + .and.(mnum.ne.5)) ind=ind+1 write (iout,'(2i5,3f10.5,5x,2e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),& @@ -4138,7 +4646,7 @@ write(iout,*) "end init MD" !el common /przechowalnia/ GGinv,gdc,Cmat,nbond !el common /przechowalnia/ nbond integer :: max_rattle = 5 - logical :: lprn = .true., lprn1 = .true., not_done + logical :: lprn = .false., lprn1 = .false., not_done real(kind=8) :: tol_rattle = 1.0d-5 integer :: nres2 nres2=2*nres @@ -4152,7 +4660,7 @@ write(iout,*) "end init MD" if (lprn) write (iout,*) "RATTLE_BROWN" nbond=nct-nnt do i=nnt,nct - if (itype(i).ne.10) nbond=nbond+1 + if (itype(i,1).ne.10) nbond=nbond+1 enddo ! Make a folded form of the Ginv-matrix ind=0 @@ -4171,7 +4679,7 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 @@ -4182,7 +4690,7 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ii=ii+1 do j=1,3 ind=ind+1 @@ -4196,7 +4704,7 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then jj=jj+1 do l=1,3 ind1=ind1+1 @@ -4228,7 +4736,7 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind1=ind1+1 do j=1,3 dC_uncor(j,ind1)=dC(j,i+nres) @@ -4244,7 +4752,7 @@ write(iout,*) "end init MD" enddo enddo do k=nnt,nct - if (itype(k).ne.10) then + if (itype(k,1).ne.10) then ind=ind+1 do j=1,3 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres) @@ -4259,9 +4767,9 @@ write(iout,*) "end init MD" x(ind)=vbld(i+1)**2-vbl**2 enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 - x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2 + x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 endif enddo if (lprn) then @@ -4278,7 +4786,7 @@ write(iout,*) "end init MD" i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t(j,i+nres),j=1,3),& @@ -4332,7 +4840,7 @@ write(iout,*) "end init MD" enddo enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 do j=1,3 xx=0.0d0 @@ -4358,9 +4866,9 @@ write(iout,*) "end init MD" i,(dC(j,i),j=1,3),x(ind),xx enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 - xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2 + xx=vbld(i+nres)**2-vbldsc0(1,itype(i,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & i,(dC(j,i+nres),j=1,3),x(ind),xx endif @@ -4373,7 +4881,7 @@ write(iout,*) "end init MD" i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i)) enddo do i=nnt,nct - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then ind=ind+1 write (iout,'(2i5,3f10.5,5x,e15.5)') & i+nres,ind,(d_t_new(j,i+nres),j=1,3),& @@ -4423,7 +4931,7 @@ write(iout,*) "end init MD" !el common /przechowalnia/ ginvfric logical :: lprn = .false., checkmode = .false. - integer :: i,j,ind,k,nres2,nres6 + integer :: i,j,ind,k,nres2,nres6,mnum nres2=2*nres nres6=6*nres @@ -4446,7 +4954,9 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + mnum=molnum(i) + if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& + .and.(mnum.ne.5)) then do j=1,3 d_t_work(ind+j)=d_t(j,i+nres) enddo @@ -4475,7 +4985,10 @@ write(iout,*) "end init MD" ind=ind+3 enddo do i=nnt,nct - if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + mnum=molnum(i) + if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& + .and.(mnum.ne.5)) then +! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 friction(j,i+nres)=fric_work(ind+j) enddo @@ -4562,6 +5075,7 @@ write(iout,*) "end init MD" !----------------------------------------------------------------------------- subroutine setup_fricmat +! use MPI use energy_data use control_data, only:time_Bcast use control, only:tcpu @@ -4595,26 +5109,30 @@ write(iout,*) "end init MD" logical :: lprn = .false. real(kind=8) :: dtdi !el ,gamvec(2*nres) !el real(kind=8),dimension(2*nres,2*nres) :: ginvfric,fcopy - real(kind=8),dimension(2*nres,2*nres) :: fcopy +! real(kind=8),allocatable,dimension(:,:) :: fcopy !el real(kind=8),dimension(2*nres*(2*nres+1)/2) :: Ghalf !(mmaxres2) (mmaxres2=(maxres2*(maxres2+1)/2)) !el common /syfek/ gamvec real(kind=8) :: work(8*2*nres) integer :: iwork(2*nres) !el common /przechowalnia/ ginvfric,Ghalf,fcopy - integer :: ii,iti,k,l,nzero,nres2,nres6,ierr + integer :: ii,iti,k,l,nzero,nres2,nres6,ierr,mnum + nres2=2*nres + nres6=6*nres #ifdef MPI + if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) + if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if (fg_rank.ne.king) goto 10 #endif - nres2=2*nres - nres6=6*nres +! nres2=2*nres +! nres6=6*nres if(.not.allocated(gamvec)) allocate(gamvec(nres2)) !(MAXRES2) if(.not.allocated(ginvfric)) allocate(ginvfric(nres2,nres2)) !maxres2=2*maxres -!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres + if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !maxres2=2*maxres !el allocate(fcopy(nres2,nres2)) !maxres2=2*maxres if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !maxres2=2*maxres -!el if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) + if(.not.allocated(fricmat)) allocate(fricmat(nres2,nres2)) ! Zeroing out fricmat do i=1,dimen do j=1,dimen @@ -4624,18 +5142,22 @@ write(iout,*) "end init MD" ! Load the friction coefficients corresponding to peptide groups ind1=0 do i=nnt,nct-1 + mnum=molnum(i) ind1=ind1+1 - gamvec(ind1)=gamp + gamvec(ind1)=gamp(mnum) enddo ! Load the friction coefficients corresponding to side chains m=nct-nnt ind=0 - gamsc(ntyp1)=1.0d0 + do j=1,2 + gamsc(ntyp1_molec(j),j)=1.0d0 + enddo do i=nnt,nct + mnum=molnum(i) ind=ind+1 ii = ind+m - iti=itype(i) - gamvec(ii)=gamsc(iabs(iti)) + iti=itype(i,mnum) + gamvec(ii)=gamsc(iabs(iti),mnum) enddo if (surfarea) call sdarea(gamvec) ! if (lprn) then @@ -4774,12 +5296,10 @@ write(iout,*) "end init MD" #else time00=tcpu() #endif -write(iout,*)"przed MPI_Scatterv in fricmat" ! Scatter the friction matrix call MPI_Scatterv(fricmat(1,1),nginv_counts(0),& nginv_start(0),MPI_DOUBLE_PRECISION,fcopy(1,1),& myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) -write(iout,*)"po MPI_Scatterv in fricmat" #ifdef TIMING #ifdef MPI time_scatter=time_scatter+MPI_Wtime()-time00 @@ -4789,13 +5309,11 @@ write(iout,*)"po MPI_Scatterv in fricmat" time_scatter_fmat=time_scatter_fmat+tcpu()-time00 #endif #endif -write(iout,*)"po MPI_Scatterv in fricmat" do i=1,dimen do j=1,2*my_ng_count fricmat(j,i)=fcopy(i,j) enddo enddo -write(iout,*)"po MPI_Scatterv in fricmat" ! write (iout,*) "My chunk of fricmat" ! call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy) endif @@ -4834,7 +5352,7 @@ write(iout,*)"po MPI_Scatterv in fricmat" real(kind=8) :: highb2,sig2,forcvec(6*nres),stochforcvec(6*nres) real(kind=8) :: time00 logical :: lprn = .false. - integer :: i,j,ind + integer :: i,j,ind,mnum do i=0,2*nres do j=1,3 @@ -4881,7 +5399,9 @@ write(iout,*)"po MPI_Scatterv in fricmat" do j=1,3 ff(j)=ff(j)+force(j,i) enddo - if (itype(i+1).ne.ntyp1) then +! if (itype(i+1,1).ne.ntyp1) then + mnum=molnum(i) + if (itype(i+1,mnum).ne.ntyp1_molec(mnum)) then do j=1,3 stochforc(j,i)=stochforc(j,i)+force(j,i+nres+1) ff(j)=ff(j)+force(j,i+nres+1) @@ -4892,7 +5412,10 @@ write(iout,*)"po MPI_Scatterv in fricmat" stochforc(j,0)=ff(j)+force(j,nnt+nres) enddo do i=nnt,nct - if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + mnum=molnum(i) + if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& + .and.(mnum.ne.5)) then +! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforc(j,i+nres)=force(j,i+nres) enddo @@ -4910,7 +5433,10 @@ write(iout,*)"po MPI_Scatterv in fricmat" ind=ind+3 enddo do i=nnt,nct - if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + mnum=molnum(i) + if ((itype(i,1).ne.10).and.(itype(i,mnum).ne.ntyp1_molec(mnum))& + .and.(mnum.ne.5)) then +! if ((itype(i,1).ne.10).and.(itype(i,1).ne.ntyp1)) then do j=1,3 stochforcvec(ind+j)=stochforc(j,i+nres) enddo @@ -5001,7 +5527,7 @@ write(iout,*)"po MPI_Scatterv in fricmat" real(kind=8),parameter :: twosix = 1.122462048309372981d0 logical :: lprn = .false. real(kind=8) :: probe,area,ratio - integer :: i,j,ind,iti + integer :: i,j,ind,iti,mnum ! ! determine new friction coefficients every few SD steps ! @@ -5015,12 +5541,14 @@ write(iout,*)"po MPI_Scatterv in fricmat" enddo ! Load peptide group radii do i=nnt,nct-1 - radius(i)=pstok + mnum=molnum(i) + radius(i)=pstok(mnum) enddo ! Load side chain radii do i=nnt,nct - iti=itype(i) - radius(i+nres)=restok(iti) + mnum=molnum(i) + iti=itype(i,mnum) + radius(i+nres)=restok(iti,mnum) enddo ! do i=1,2*nres ! write (iout,*) "i",i," radius",radius(i) @@ -5046,7 +5574,8 @@ write(iout,*)"po MPI_Scatterv in fricmat" ind=ind+1 gamvec(ind) = ratio * gamvec(ind) enddo - stdforcp(i)=stdfp*dsqrt(gamvec(ind)) + mnum=molnum(i) + stdforcp(i)=stdfp(mnum)*dsqrt(gamvec(ind)) if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i) endif enddo @@ -5060,7 +5589,8 @@ write(iout,*)"po MPI_Scatterv in fricmat" ind=ind+1 gamvec(ind) = ratio * gamvec(ind) enddo - stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind)) + mnum=molnum(i) + stdforcsc(i)=stdfsc(itype(i,mnum),mnum)*dsqrt(gamvec(ind)) if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i) endif enddo @@ -5165,11 +5695,11 @@ write(iout,*)"po MPI_Scatterv in fricmat" logical :: omit(maxarc) ! ! include 'sizes.i' - maxatm = 2*nres !maxres2 maxres2=2*maxres - maxlight = 8*maxatm - maxbnd = 2*maxatm - maxang = 3*maxatm - maxtors = 4*maxatm +! maxatm = 2*nres !maxres2 maxres2=2*maxres +! maxlight = 8*maxatm +! maxbnd = 2*maxatm +! maxang = 3*maxatm +! maxtors = 4*maxatm ! ! zero out the surface area for the sphere of interest ! @@ -5610,7 +6140,7 @@ write(iout,*)"po MPI_Scatterv in fricmat" allocate(flag_stoch(0:maxflag_stoch)) !(0:maxflag_stoch) #endif -!el if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) + if(.not.allocated(fcopy)) allocate(fcopy(nres2,nres2)) !---------------------- ! commom.hairpin in CSA module !---------------------- @@ -5626,9 +6156,14 @@ write(iout,*)"po MPI_Scatterv in fricmat" ! common /lagrange/ allocate(d_t(3,0:nres2),d_a(3,0:nres2),d_t_old(3,0:nres2)) !(3,0:MAXRES2) allocate(d_a_work(nres6)) !(6*MAXRES) +#ifdef FIVEDIAG + allocate(DM(nres2),DU1(nres2),DU2(nres2)) + allocate(DMorig(nres2),DU1orig(nres2),DU2orig(nres2)) +#else allocate(Gmat(nres2,nres2),A(nres2,nres2)) if(.not.allocated(Ginv)) allocate(Ginv(nres2,nres2)) !in control: ergastulum allocate(Gsqrp(nres2,nres2),Gsqrm(nres2,nres2),Gvec(nres2,nres2)) !(maxres2,maxres2) +#endif allocate(Geigen(nres2)) !(maxres2) if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2) ! common /inertia/ in io_conf: parmread @@ -5645,6 +6180,14 @@ write(iout,*)"po MPI_Scatterv in fricmat" if(.not.allocated(nginv_start)) allocate(nginv_start(0:nfgtasks)) !(0:MaxProcs) !---------------------- ! common.muca in read_muca +! common /double_muca/ +! real(kind=8) :: elow,ehigh,factor,hbin,factor_min +! real(kind=8),dimension(:),allocatable :: emuca,nemuca,& +! nemuca2,hist !(4*maxres) +! common /integer_muca/ +! integer :: nmuca,imtime,muca_smooth +! common /mucarem/ +! real(kind=8),dimension(:),allocatable :: elowi,ehighi !(maxprocs) !---------------------- ! common.MD ! common /mdgrad/ in module.energy