X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FREMD.f90;h=e2db47896619eb6a519ebd30fb2bccab65fb94e3;hb=705644e0cbb7678faefd6fe1bc436159d38ad85d;hp=1e294adfc1047114a57826de290044372ce756cf;hpb=8ca97b16fe25b7053f258263899ba030572cc58f;p=unres4.git diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90 index 1e294ad..e2db478 100644 --- a/source/unres/REMD.f90 +++ b/source/unres/REMD.f90 @@ -47,7 +47,7 @@ #ifdef CHECK5DSOL real(kind=8) :: rscheck(dimen),rsold(dimen) #endif - logical :: lprn = .true. + logical :: lprn = .false. !el common /cipiszcze/ itime itime = itt_comm @@ -62,14 +62,16 @@ enddo write (iout,*) "Potential forces sidechain" do i=nnt,nct - if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) & + mnum=molnum(i) + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.eq.5)& write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3) enddo endif do j=1,3 ind=1 do i=nnt,nct - if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then + mnum=molnum(i) + if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.eq.5)then rs(ind)=-gcart(j,i)-gxcart(j,i) ind=ind+1 else @@ -192,7 +194,8 @@ if (lprn) write (iout,*) "Potential forces sidechain" do i=nnt,nct 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 if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') & i,(-gxcart(j,i),j=1,3) do j=1,3 @@ -217,7 +220,8 @@ do i=nnt,nct 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)) 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_a(j,i+nres)=d_a_work(ind) @@ -296,10 +300,10 @@ !#endif ! include 'COMMON.IOUNITS' ! include 'COMMON.TIME1' - logical :: lprn = .true. + logical :: lprn = .false. logical :: osob #ifndef FIVEDIAG - real(kind=8),dimension(2*nres,2*nres) :: Bmat,matmult + real(kind=8),allocatable,dimension(:,:) :: Bmat,matmult #endif real(kind=8) :: dtdi real(kind=8),dimension(2*nres) :: massvec,sqreig !(maxres2) maxres2=2*maxres @@ -311,18 +315,22 @@ real(kind=8),dimension(8*6*nres) :: work !(8*maxres6) integer,dimension(6*nres) :: iwork !(maxres6) maxres6=6*maxres !el common /przechowalnia/ Gcopy,Ghalf - real(kind=8) :: coeff + real(kind=8) :: coeff,mscab integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark real(kind=8) :: ip4 integer :: iz,mnum + print *,"just entered" relfeh=1.0d-14 nres2=2*nres - + print *,"FIVEDIAG" write (iout,*) "before FIVEDIAG" #ifndef FIVEDIAG write (iout,*) "ALLOCATE" + print *,"ALLOCATE" if(.not.allocated(Gcopy)) allocate(Gcopy(nres2,nres2)) !(maxres2,maxres2) if(.not.allocated(Ghalf)) allocate(Ghalf(nres2*(nres2+1)/2)) !mmaxres2=(maxres2*(maxres+1)/2) + if(.not.allocated(Bmat)) allocate(Bmat(nres2,nres2)) + if(.not.allocated(matmult)) allocate(matmult(nres2,nres2)) #endif ! ! Set up the matrix of the (dC,dX)-->(C,X) transformation (A), the @@ -335,6 +343,7 @@ dimen3=dimen*3 write (iout,*) "nnt",nnt," nct",nct," nside",nside #ifdef FIVEDIAG +#ifdef CRYST_BOND ip4=ip/4 DM(1)=mp/4+ip4 if (iabs(itype(nnt)).eq.10) then @@ -374,7 +383,9 @@ ind=1 do i=nnt,nct - if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then + mnum=molnum(i) + if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 & + .and.mnum.eq.5) then DU1(ind)=-isc(iabs(itype(i,1))) DU1(ind+1)=0.0d0 ind=ind+2 @@ -400,6 +411,75 @@ ind=ind+1 endif enddo +#else + ip4=ip(1)/4 + DM(1)=mp(1)/4+ip4 + if (iabs(itype(nnt,1)).eq.10) then + DM(1)=DM(1)+msc(10,1) + ind=2 + nind=1 + else + DM(1)=DM(1)+isc(iabs(itype(nnt,1)),1) + DM(2)=msc(iabs(itype(nnt,1)),1)+isc(iabs(itype(nnt,1)),1) + ind=3 + nind=2 + endif + do i=nnt+1,nct-1 +! if (iabs(itype(i,1)).eq.ntyp1) cycle + DM(ind)=2*ip4+mp(1)/2 + if (iabs(itype(i,1)).eq.10 .or. iabs(itype(i,1)).eq.ntyp1) then + if (iabs(itype(i,1)).eq.10) DM(ind)=DM(ind)+msc(10,1) + ind=ind+1 + else + DM(ind)=DM(ind)+isc(iabs(itype(i,1)),1) + DM(ind+1)=msc(iabs(itype(i,1)),1)+isc(iabs(itype(i,1)),1) + ind=ind+2 + endif + enddo + if (nct.gt.nnt) then + DM(ind)=ip4+mp(1)/4 + if (iabs(itype(nct,1)).eq.10) then + DM(ind)=DM(ind)+msc(10,1) + nind=ind + else + DM(ind)=DM(ind)+isc(iabs(itype(nct,1)),1) + DM(ind+1)=msc(iabs(itype(nct,1)),1)+isc(iabs(itype(nct,1)),1) + nind=ind+1 + endif + endif + + + ind=1 + do i=nnt,nct + mnum=molnum(i) + if (iabs(itype(i,1)).ne.10 .and.iabs(itype(i,mnum)).ne.ntyp1 & + .and.mnum.eq.5) then + DU1(ind)=-isc(iabs(itype(i,1)),1) + DU1(ind+1)=0.0d0 + ind=ind+2 + else + DU1(ind)=mp(1)/4-ip4 + ind=ind+1 + endif + enddo + + ind=1 + do i=nnt,nct-1 + mnum=molnum(i) +! if (iabs(itype(i,1)).eq.ntyp1) cycle + write (iout,*) "i",i," itype",itype(i,1),ntyp1 + if (iabs(itype(i,1)).ne.10 .and. & + iabs(itype(i,mnum)).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then + DU2(ind)=mp(1)/4-ip4 + DU2(ind+1)=0.0d0 + ind=ind+2 + else + DU2(ind)=0.0d0 + DU2(ind+1)=0.0d0 + ind=ind+1 + endif + enddo +#endif DMorig=DM DU1orig=DU1 DU2orig=DU2 @@ -513,6 +593,8 @@ ind=ind+1 ind1=ind1+1 coeff=0.25d0*IP(mnum) + if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum) + print *,i,coeff,mp(mnum) massvec(ind1)=mp(mnum) Gmat(ind,ind)=coeff print *,"i",mp(mnum) @@ -539,7 +621,13 @@ ii = ind+m mnum=molnum(i) iti=itype(i,mnum) - massvec(ii)=msc(iabs(iti),mnum) + if (mnum.eq.5) then + mscab=0.0 + else + mscab=msc(iabs(iti),mnum) + endif + massvec(ii)=mscab + if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then ind1=ind1+1 ii1= ind1+m1 @@ -627,15 +715,20 @@ call matout(dimen,dimen,nres2,nres2,Bmat) endif Gcopy=0.0d0 + write(iout,*) "before Gcopy",dimen,nres*2 +#ifdef TESTFIVEDIAG do i=1,dimen do j=1,dimen +! write (iout,*) "myij",i,j do k=1,dimen do l=1,dimen +! write(iout,*) "i,j,k,l",i,j,k,l Gcopy(i,j)=Gcopy(i,j)+Bmat(k,i)*Gmat(k,l)*Bmat(l,j) enddo enddo enddo enddo +#endif if (lprn) then write (iout,'(//a)') "Gmat-transformed" call matout(dimen,dimen,nres2,nres2,Gcopy)