X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=1d44fc1b2545f9da62f5291090ff5a16616060b8;hb=a9930dd81fb12f9ca37898d19e63d8064bf91e03;hp=3d3dc39ec738a7635562d24521a674cb9963ae72;hpb=4367d241fbb2bc284580092d2d177b7c79ac3a42;p=unres4.git diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 3d3dc39..1d44fc1 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -141,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) !----------------------------------------------------------------------------- ! ! @@ -184,9 +184,12 @@ integer :: rstcount !ilen, !el external ilen !el real(kind=8),dimension(6*nres) :: stochforcvec !(MAXRES6) maxres6=6*maxres - 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(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 @@ -202,6 +205,9 @@ logical :: osob nres2=2*nres nres6=6*nres + allocate(Bmat(6*nres,2*nres),GBmat(6*nres,2*nres),Tmat(6*nres,2*nres)) + allocate(Cmat_(2*nres,2*nres),Cinv(2*nres,2*nres)) + allocate(Pmat(6*nres,6*nres)) if (.not.allocated(stochforcvec)) allocate(stochforcvec(nres6)) !(MAXRES6) maxres6=6*maxres @@ -652,9 +658,9 @@ ! 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 @@ -672,12 +678,14 @@ 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 @@ -690,8 +698,16 @@ incr(j)=d_t(j,0) enddo do i=nnt,nct - iti=iabs(itype(i,1)) - if (itype(i,1).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 @@ -702,7 +718,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) @@ -713,13 +729,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 @@ -727,20 +744,23 @@ ! The rotational part of the side chain virtual bond KEr_sc=0.0D0 do i=nnt,nct - iti=iabs(itype(i,1)) - if (itype(i,1).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 @@ -796,8 +816,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') @@ -816,10 +838,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 @@ -863,7 +889,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 @@ -956,7 +984,10 @@ #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 + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -969,7 +1000,8 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -990,6 +1022,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) @@ -1767,7 +1800,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 @@ -1778,7 +1811,10 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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 @@ -1804,7 +1840,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" @@ -1830,7 +1866,10 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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 @@ -1866,7 +1905,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 @@ -1877,7 +1916,11 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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 @@ -1960,7 +2003,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 @@ -1984,7 +2027,11 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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 @@ -2022,7 +2069,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. @@ -2052,7 +2099,11 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) & @@ -2084,7 +2135,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) @@ -2122,7 +2173,11 @@ enddo endif do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -2179,7 +2234,8 @@ enddo endif ! Side chains - if (itype(i,1).ne.10 .and. itype(i,1).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)) @@ -2212,7 +2268,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)) @@ -2228,7 +2284,11 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -2283,7 +2343,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 @@ -2293,11 +2353,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,1))) & - *dsqrt(gamsc(iabs(itype(i,1)))) + 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 @@ -2649,7 +2711,7 @@ real (kind=8),allocatable,dimension(:,:) :: matold #endif #endif - integer :: i,j,ii,k,ind,mark,imark + 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 @@ -2873,7 +2935,13 @@ Ek1=0.0d0 ii=0 do i=nnt,nct - if (itype(i,1).eq.10) then +! 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 @@ -2881,18 +2949,21 @@ 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*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& - 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 + 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) ii=ii+3 - write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1)) + 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 (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1)))*(d_t_work(ii)-d_t_work(ii-3))**2 - Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)))*d_t_work(ii)**2 + 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 @@ -2916,7 +2987,9 @@ d_t(k,j)=d_t_work(ind) ind=ind+1 enddo - if (itype(j).ne.10 .and. itype(j).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 k=1,3 d_t(k,j+nres)=d_t_work(ind) ind=ind+1 @@ -2936,7 +3009,10 @@ d_t(j,0)=d_t(j,nnt) enddo do i=nnt,nct - if (itype(i,1).eq.10) then +! 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 @@ -2984,7 +3060,10 @@ enddo enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -3225,7 +3304,10 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -3273,7 +3355,10 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).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) @@ -3339,7 +3424,10 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -3552,7 +3640,10 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -3618,7 +3709,10 @@ ind=ind+3 enddo do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).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) @@ -3654,13 +3748,13 @@ 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 @@ -3672,84 +3766,99 @@ 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,1)) - 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,1)) + 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,1).ne.10 .and. itype(i,1).ne.ntyp1) then - iti=iabs(itype(i,1)) + 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 @@ -3819,7 +3928,11 @@ enddo enddo do i=nnt,nct - if(itype(i,1).ne.10 .and. itype(i,1).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 @@ -3849,9 +3962,9 @@ ! 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 @@ -3860,6 +3973,8 @@ 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 @@ -3871,7 +3986,7 @@ 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) @@ -3879,19 +3994,26 @@ 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,1)) + 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,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 v(j)=incr(j)+d_t(j,inres) enddo @@ -3904,16 +4026,17 @@ ! 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,1).ne.10 .and. itype(i,1).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 @@ -3938,7 +4061,7 @@ ! 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 @@ -3946,15 +4069,22 @@ 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,1))) + if (mnum.ne.5) then + amas=msc(iabs(itype(i,mnum)),mnum) + else + amas=0.0d0 + endif summas=summas+amas - if (itype(i,1).ne.10 .and. itype(i,1).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 @@ -4023,7 +4153,9 @@ if (lprn) write (iout,*) "RATTLE1" nbond=nct-nnt do i=nnt,nct - if (itype(i,1).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 @@ -4042,7 +4174,9 @@ enddo enddo do k=nnt,nct - if (itype(k,1).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 @@ -4053,7 +4187,9 @@ enddo enddo do i=nnt,nct - if (itype(i,1).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 @@ -4149,7 +4285,10 @@ 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,1).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),& @@ -4229,7 +4368,9 @@ i,(dC(j,i),j=1,3),x(ind),xx enddo do i=nnt,nct - if (itype(i,1).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,1))**2 write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') & @@ -4319,7 +4460,9 @@ enddo enddo do k=nnt,nct - if (itype(k,1).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) @@ -4348,7 +4491,9 @@ enddo enddo do i=nnt,nct - if (itype(i,1).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 @@ -4365,7 +4510,9 @@ x(ind)=scalar(d_t(1,i),dC(1,i)) enddo do i=nnt,nct - if (itype(i,1).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 @@ -4379,7 +4526,9 @@ i,ind,(d_t(j,i),j=1,3),x(ind) enddo do i=nnt,nct - if (itype(i,1).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) @@ -4414,7 +4563,9 @@ enddo enddo do i=nnt,nct - if (itype(i,1).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 @@ -4435,7 +4586,9 @@ 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,1).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),& @@ -4772,7 +4925,7 @@ !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 @@ -4795,7 +4948,9 @@ ind=ind+3 enddo do i=nnt,nct - if ((itype(i,1).ne.10).and.(itype(i,1).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 @@ -4824,7 +4979,10 @@ ind=ind+3 enddo do i=nnt,nct - if ((itype(i,1).ne.10).and.(itype(i,1).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 @@ -4945,26 +5103,30 @@ 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 @@ -4974,18 +5136,22 @@ ! 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,1) - gamvec(ii)=gamsc(iabs(iti)) + iti=itype(i,mnum) + gamvec(ii)=gamsc(iabs(iti),mnum) enddo if (surfarea) call sdarea(gamvec) ! if (lprn) then @@ -5180,7 +5346,7 @@ 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 @@ -5227,7 +5393,9 @@ do j=1,3 ff(j)=ff(j)+force(j,i) enddo - if (itype(i+1,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) @@ -5238,7 +5406,10 @@ stochforc(j,0)=ff(j)+force(j,nnt+nres) enddo do i=nnt,nct - if ((itype(i,1).ne.10).and.(itype(i,1).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 @@ -5256,7 +5427,10 @@ ind=ind+3 enddo do i=nnt,nct - if ((itype(i,1).ne.10).and.(itype(i,1).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 @@ -5347,7 +5521,7 @@ 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 ! @@ -5361,12 +5535,14 @@ 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,1) - 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) @@ -5392,7 +5568,8 @@ 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 @@ -5406,7 +5583,8 @@ ind=ind+1 gamvec(ind) = ratio * gamvec(ind) enddo - stdforcsc(i)=stdfsc(itype(i,1))*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 @@ -5956,7 +6134,7 @@ 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 !----------------------