+++ /dev/null
-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
projects / unres.git / blobdiff