X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fmd-diff%2Fmd%2Fbrown_step.f;fp=source%2Funres%2Fsrc_MD%2Fmd-diff%2Fmd%2Fbrown_step.f;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=050aac869c4290b73bfd7a1980e7c397cc119c80;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/md-diff/md/brown_step.f b/source/unres/src_MD/md-diff/md/brown_step.f deleted file mode 100644 index 050aac8..0000000 --- a/source/unres/src_MD/md-diff/md/brown_step.f +++ /dev/null @@ -1,377 +0,0 @@ -c------------------------------------------------------------------------------- - subroutine brown_step(itime) -c------------------------------------------------ -c Perform a single Euler integration step of Brownian dynamics -c------------------------------------------------ - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CONTROL' - include 'COMMON.VAR' - include 'COMMON.MD' - include 'COMMON.LANGEVIN' - include 'COMMON.CHAIN' - include 'COMMON.DERIV' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.INTERACT' - include 'COMMON.IOUNITS' - include 'COMMON.NAMES' - include 'COMMON.TIME1' - double precision energia(0:n_ene),zapas(MAXRES6) - integer ilen,rstcount - external ilen - double precision stochforcvec(MAXRES6) - double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2), - & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2), - & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6), - & ppvec(maxres2) - common /stochcalc/ stochforcvec - common /gucio/ cm, energia - integer itime - logical lprn /.false./,lprn1 /.false./ - integer maxiter /5/ - double precision difftol /1.0d-5/ - nbond=nct-nnt - do i=nnt,nct - if (itype(i).ne.10) nbond=nbond+1 - enddo -c - if (lprn1) then - write (iout,*) "Generalized inverse of fricmat" - call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat) - endif - do i=1,dimen - do j=1,nbond - Bmat(i,j)=0.0d0 - enddo - enddo - ind=3 - ind1=0 - do i=nnt,nct-1 - ind1=ind1+1 - do j=1,3 - Bmat(ind+j,ind1)=dC_norm(j,i) - enddo - ind=ind+3 - enddo - do i=nnt,nct - if (itype(i).ne.10) then - ind1=ind1+1 - do j=1,3 - Bmat(ind+j,ind1)=dC_norm(j,i+nres) - enddo - ind=ind+3 - endif - enddo - if (lprn1) then - write (iout,*) "Matrix Bmat" - call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat) - endif - do i=1,dimen - do j=1,nbond - GBmat(i,j)=0.0d0 - do k=1,dimen - GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j) - enddo - enddo - enddo - if (lprn1) then - write (iout,*) "Matrix GBmat" - call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat) - endif - do i=1,nbond - do j=1,nbond - Cmat(i,j)=0.0d0 - do k=1,dimen - Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j) - enddo - enddo - enddo - if (lprn1) then - write (iout,*) "Matrix Cmat" - call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat) - endif - call matinvert(nbond,MAXRES2,Cmat,Cinv) - if (lprn1) then - write (iout,*) "Matrix Cinv" - call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv) - endif - do i=1,dimen - do j=1,nbond - Tmat(i,j)=0.0d0 - do k=1,nbond - Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j) - enddo - enddo - enddo - if (lprn1) then - write (iout,*) "Matrix Tmat" - call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat) - endif - do i=1,dimen - do j=1,dimen - if (i.eq.j) then - Pmat(i,j)=1.0d0 - else - Pmat(i,j)=0.0d0 - endif - do k=1,nbond - Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k) - enddo - enddo - enddo - if (lprn1) then - write (iout,*) "Matrix Pmat" - call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat) - endif - do i=1,dimen - Td(i)=0.0d0 - ind=0 - do k=nnt,nct-1 - ind=ind+1 - Td(i)=Td(i)+vbl*Tmat(i,ind) - enddo - do k=nnt,nct - if (itype(k).ne.10) then - ind=ind+1 - Td(i)=Td(i)+vbldsc0(itype(k))*Tmat(i,ind) - endif - enddo - enddo - if (lprn1) then - write (iout,*) "Vector Td" - do i=1,dimen - write (iout,'(i5,f10.5)') i,Td(i) - enddo - endif - call stochastic_force(stochforcvec) - if (lprn) then - write (iout,*) "stochforcvec" - do i=1,dimen - write (iout,*) i,stochforcvec(i) - enddo - endif - do j=1,3 - zapas(j)=-gcart(j,0)+stochforcvec(j) - d_t_work(j)=d_t(j,0) - dC_work(j)=dC_old(j,0) - enddo - ind=3 - do i=nnt,nct-1 - do j=1,3 - ind=ind+1 - zapas(ind)=-gcart(j,i)+stochforcvec(ind) - dC_work(ind)=dC_old(j,i) - enddo - enddo - do i=nnt,nct - if (itype(i).ne.10) then - do j=1,3 - ind=ind+1 - zapas(ind)=-gxcart(j,i)+stochforcvec(ind) - dC_work(ind)=dC_old(j,i+nres) - enddo - endif - enddo - - if (lprn) then - write (iout,*) "Initial d_t_work" - do i=1,dimen - write (iout,*) i,d_t_work(i) - enddo - endif - - do i=1,dimen - d_t_work(i)=0.0d0 - do j=1,dimen - d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j) - enddo - enddo - - do i=1,dimen - zapas(i)=Td(i) - do j=1,dimen - zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time) - enddo - enddo - if (lprn1) then - write (iout,*) "Final d_t_work and zapas" - do i=1,dimen - write (iout,*) i,d_t_work(i),zapas(i) - enddo - endif - - do j=1,3 - d_t(j,0)=d_t_work(j) - dc(j,0)=zapas(j) - dc_work(j)=dc(j,0) - enddo - ind=3 - do i=nnt,nct-1 - do j=1,3 - d_t(j,i)=d_t_work(i) - dc(j,i)=zapas(ind+j) - dc_work(ind+j)=dc(j,i) - enddo - ind=ind+3 - enddo - do i=nnt,nct - do j=1,3 - d_t(j,i+nres)=d_t_work(ind+j) - dc(j,i+nres)=zapas(ind+j) - dc_work(ind+j)=dc(j,i+nres) - enddo - ind=ind+3 - enddo - if (lprn) then - call chainbuild_cart - write (iout,*) "Before correction for rotational lengthening" - write (iout,*) "New coordinates", - & " and differences between actual and standard bond lengths" - ind=0 - do i=nnt,nct-1 - ind=ind+1 - xx=vbld(i+1)-vbl - write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') - & i,(dC(j,i),j=1,3),xx - enddo - do i=nnt,nct - if (itype(i).ne.10) then - ind=ind+1 - xx=vbld(i+nres)-vbldsc0(itype(i)) - write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') - & i,(dC(j,i+nres),j=1,3),xx - endif - enddo - endif -c Second correction (rotational lengthening) -c do iter=1,maxiter - diffmax=0.0d0 - ind=0 - do i=nnt,nct-1 - ind=ind+1 - blen2 = scalar(dc(1,i),dc(1,i)) - ppvec(ind)=2*vbl**2-blen2 - diffbond=dabs(vbl-dsqrt(blen2)) - if (diffbond.gt.diffmax) diffmax=diffbond - if (ppvec(ind).gt.0.0d0) then - ppvec(ind)=dsqrt(ppvec(ind)) - else - ppvec(ind)=0.0d0 - endif - if (lprn) then - write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) - endif - enddo - do i=nnt,nct - if (itype(i).ne.10) then - ind=ind+1 - blen2 = scalar(dc(1,i+nres),dc(1,i+nres)) - ppvec(ind)=2*vbldsc0(itype(i))**2-blen2 - diffbond=dabs(vbldsc0(itype(i))-dsqrt(blen2)) - if (diffbond.gt.diffmax) diffmax=diffbond - if (ppvec(ind).gt.0.0d0) then - ppvec(ind)=dsqrt(ppvec(ind)) - else - ppvec(ind)=0.0d0 - endif - if (lprn) then - write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind) - endif - endif - enddo - if (lprn) write (iout,*) "iter",iter," diffmax",diffmax - if (diffmax.lt.difftol) goto 10 - do i=1,dimen - Td(i)=0.0d0 - do j=1,nbond - Td(i)=Td(i)+ppvec(j)*Tmat(i,j) - enddo - enddo - do i=1,dimen - zapas(i)=Td(i) - do j=1,dimen - zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j) - enddo - enddo - do j=1,3 - dc(j,0)=zapas(j) - dc_work(j)=zapas(j) - enddo - ind=3 - do i=nnt,nct-1 - do j=1,3 - dc(j,i)=zapas(ind+j) - dc_work(ind+j)=zapas(ind+j) - enddo - ind=ind+3 - enddo - do i=nnt,nct - if (itype(i).ne.10) then - do j=1,3 - dc(j,i+nres)=zapas(ind+j) - dc_work(ind+j)=zapas(ind+j) - enddo - ind=ind+3 - endif - enddo -c Building the chain from the newly calculated coordinates - call chainbuild_cart - if (large.and. mod(itime,ntwe).eq.0) then - write (iout,*) "Cartesian and internal coordinates: step 1" - call cartprint - call intout - write (iout,'(a)') "Potential forces" - do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3), - & (-gxcart(j,i),j=1,3) - enddo - write (iout,'(a)') "Stochastic forces" - do i=0,nres - write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3), - & (stochforc(j,i+nres),j=1,3) - enddo - write (iout,'(a)') "Velocities" - 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 - if (lprn) then - write (iout,*) "After correction for rotational lengthening" - write (iout,*) "New coordinates", - & " and differences between actual and standard bond lengths" - ind=0 - do i=nnt,nct-1 - ind=ind+1 - xx=vbld(i+1)-vbl - write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') - & i,(dC(j,i),j=1,3),xx - enddo - do i=nnt,nct - if (itype(i).ne.10) then - ind=ind+1 - xx=vbld(i+nres)-vbldsc0(itype(i)) - write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') - & i,(dC(j,i+nres),j=1,3),xx - endif - enddo - endif -c ENDDO -c write (iout,*) "Too many attempts at correcting the bonds" -c stop - 10 continue - tt0 = tcpu() -c Calculate energy and forces - call zerograd - call etotal(energia) - potE=energia(0)-energia(20) - call cartgrad - totT=totT+d_time -c Calculate the kinetic and total energy and the kinetic temperature - call kinetic(EK) - t_enegrad=t_enegrad+tcpu()-tt0 - totE=EK+potE - kinetic_T=2.0d0/(dimen*Rb)*EK - return - end