X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fstochfric.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fstochfric.F;h=13d02fb43d59eaf94cdc84a0290b5a4e290d1d8a;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/stochfric.F b/source/unres/src_MD-M-newcorr/stochfric.F new file mode 100644 index 0000000..13d02fb --- /dev/null +++ b/source/unres/src_MD-M-newcorr/stochfric.F @@ -0,0 +1,627 @@ + subroutine friction_force + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.MD' +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.IOUNITS' + double precision gamvec(MAXRES6) + common /syfek/ gamvec + double precision vv(3),vvtot(3,maxres),v_work(MAXRES6), + & ginvfric(maxres2,maxres2) + common /przechowalnia/ ginvfric + + logical lprn /.false./, checkmode /.false./ + + do i=0,MAXRES2 + do j=1,3 + friction(j,i)=0.0d0 + enddo + enddo + + do j=1,3 + d_t_work(j)=d_t(j,0) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + d_t_work(ind+j)=d_t(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + do j=1,3 + d_t_work(ind+j)=d_t(j,i+nres) + enddo + ind=ind+3 + endif + enddo + + call fricmat_mult(d_t_work,fric_work) + + if (.not.checkmode) return + + if (lprn) then + write (iout,*) "d_t_work and fric_work" + do i=1,3*dimen + write (iout,'(i3,2e15.5)') i,d_t_work(i),fric_work(i) + enddo + endif + do j=1,3 + friction(j,0)=fric_work(j) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + friction(j,i)=fric_work(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + do j=1,3 + friction(j,i+nres)=fric_work(ind+j) + enddo + ind=ind+3 + endif + enddo + if (lprn) then + write(iout,*) "Friction backbone" + do i=0,nct-1 + write(iout,'(i5,3e15.5,5x,3e15.5)') + & i,(friction(j,i),j=1,3),(d_t(j,i),j=1,3) + enddo + write(iout,*) "Friction side chain" + do i=nnt,nct + write(iout,'(i5,3e15.5,5x,3e15.5)') + & i,(friction(j,i+nres),j=1,3),(d_t(j,i+nres),j=1,3) + enddo + endif + if (lprn) then + do j=1,3 + vv(j)=d_t(j,0) + enddo + do i=nnt,nct + do j=1,3 + vvtot(j,i)=vv(j)+0.5d0*d_t(j,i) + vvtot(j,i+nres)=vv(j)+d_t(j,i+nres) + vv(j)=vv(j)+d_t(j,i) + enddo + enddo + write (iout,*) "vvtot backbone and sidechain" + do i=nnt,nct + write (iout,'(i5,3e15.5,5x,3e15.5)') i,(vvtot(j,i),j=1,3), + & (vvtot(j,i+nres),j=1,3) + enddo + ind=0 + do i=nnt,nct-1 + do j=1,3 + v_work(ind+j)=vvtot(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + do j=1,3 + v_work(ind+j)=vvtot(j,i+nres) + enddo + ind=ind+3 + enddo + write (iout,*) "v_work gamvec and site-based friction forces" + do i=1,dimen1 + write (iout,'(i5,3e15.5)') i,v_work(i),gamvec(i), + & gamvec(i)*v_work(i) + enddo +c do i=1,dimen +c fric_work1(i)=0.0d0 +c do j=1,dimen1 +c fric_work1(i)=fric_work1(i)-A(j,i)*gamvec(j)*v_work(j) +c enddo +c enddo +c write (iout,*) "fric_work and fric_work1" +c do i=1,dimen +c write (iout,'(i5,2e15.5)') i,fric_work(i),fric_work1(i) +c enddo + do i=1,dimen + do j=1,dimen + ginvfric(i,j)=0.0d0 + do k=1,dimen + ginvfric(i,j)=ginvfric(i,j)+ginv(i,k)*fricmat(k,j) + enddo + enddo + enddo + write (iout,*) "ginvfric" + do i=1,dimen + write (iout,'(i5,100f8.3)') i,(ginvfric(i,j),j=1,dimen) + enddo + write (iout,*) "symmetry check" + do i=1,dimen + do j=1,i-1 + write (iout,*) i,j,ginvfric(i,j)-ginvfric(j,i) + enddo + enddo + endif + return + end +c----------------------------------------------------- + subroutine stochastic_force(stochforcvec) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.MD' + include 'COMMON.TIME1' +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.IOUNITS' + + double precision x,sig,lowb,highb, + & ff(3),force(3,0:MAXRES2),zeta2,lowb2, + & highb2,sig2,forcvec(MAXRES6),stochforcvec(MAXRES6) + logical lprn /.false./ + do i=0,MAXRES2 + do j=1,3 + stochforc(j,i)=0.0d0 + enddo + enddo + x=0.0d0 + +#ifdef MPI + time00=MPI_Wtime() +#else + time00=tcpu() +#endif +c Compute the stochastic forces acting on bodies. Store in force. + do i=nnt,nct-1 + sig=stdforcp(i) + lowb=-5*sig + highb=5*sig + do j=1,3 + force(j,i)=anorm_distr(x,sig,lowb,highb) + enddo + enddo + do i=nnt,nct + sig2=stdforcsc(i) + lowb2=-5*sig2 + highb2=5*sig2 + do j=1,3 + force(j,i+nres)=anorm_distr(x,sig2,lowb2,highb2) + enddo + enddo +#ifdef MPI + time_fsample=time_fsample+MPI_Wtime()-time00 +#else + time_fsample=time_fsample+tcpu()-time00 +#endif +c Compute the stochastic forces acting on virtual-bond vectors. + do j=1,3 + ff(j)=0.0d0 + enddo + do i=nct-1,nnt,-1 + do j=1,3 + stochforc(j,i)=ff(j)+0.5d0*force(j,i) + enddo + do j=1,3 + ff(j)=ff(j)+force(j,i) + enddo + if (itype(i+1).ne.ntyp1) 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) + enddo + endif + enddo + do j=1,3 + 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 + do j=1,3 + stochforc(j,i+nres)=force(j,i+nres) + enddo + endif + enddo + + do j=1,3 + stochforcvec(j)=stochforc(j,0) + enddo + ind=3 + do i=nnt,nct-1 + do j=1,3 + stochforcvec(ind+j)=stochforc(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then + do j=1,3 + stochforcvec(ind+j)=stochforc(j,i+nres) + enddo + ind=ind+3 + endif + enddo + if (lprn) then + write (iout,*) "stochforcvec" + do i=1,3*dimen + write(iout,'(i5,e15.5)') i,stochforcvec(i) + enddo + write(iout,*) "Stochastic forces backbone" + do i=0,nct-1 + write(iout,'(i5,3e15.5)') i,(stochforc(j,i),j=1,3) + enddo + write(iout,*) "Stochastic forces side chain" + do i=nnt,nct + write(iout,'(i5,3e15.5)') + & i,(stochforc(j,i+nres),j=1,3) + enddo + endif + + if (lprn) then + + ind=0 + do i=nnt,nct-1 + write (iout,*) i,ind + do j=1,3 + forcvec(ind+j)=force(j,i) + enddo + ind=ind+3 + enddo + do i=nnt,nct + write (iout,*) i,ind + do j=1,3 + forcvec(j+ind)=force(j,i+nres) + enddo + ind=ind+3 + enddo + + write (iout,*) "forcvec" + ind=0 + do i=nnt,nct-1 + do j=1,3 + write (iout,'(2i3,2f10.5)') i,j,force(j,i), + & forcvec(ind+j) + enddo + ind=ind+3 + enddo + do i=nnt,nct + do j=1,3 + write (iout,'(2i3,2f10.5)') i,j,force(j,i+nres), + & forcvec(ind+j) + enddo + ind=ind+3 + enddo + + endif + + return + end +c------------------------------------------------------------------ + subroutine setup_fricmat + implicit real*8 (a-h,o-z) +#ifdef MPI + include 'mpif.h' +#endif + include 'DIMENSIONS' + include 'COMMON.VAR' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.MD' + include 'COMMON.SETUP' + include 'COMMON.TIME1' +c integer licznik /0/ +c save licznik +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.IOUNITS' + integer IERROR + integer i,j,ind,ind1,m + logical lprn /.false./ + double precision dtdi,gamvec(MAXRES2), + & ginvfric(maxres2,maxres2),Ghalf(mmaxres2),fcopy(maxres2,maxres2) + common /syfek/ gamvec + double precision work(8*maxres2) + integer iwork(maxres2) + common /przechowalnia/ ginvfric,Ghalf,fcopy +#ifdef MPI + if (fg_rank.ne.king) goto 10 +#endif +c Zeroing out fricmat + do i=1,dimen + do j=1,dimen + fricmat(i,j)=0.0d0 + enddo + enddo +c Load the friction coefficients corresponding to peptide groups + ind1=0 + do i=nnt,nct-1 + ind1=ind1+1 + gamvec(ind1)=gamp + enddo +c Load the friction coefficients corresponding to side chains + m=nct-nnt + ind=0 + gamsc(ntyp1)=1.0d0 + do i=nnt,nct + ind=ind+1 + ii = ind+m + iti=itype(i) + gamvec(ii)=gamsc(iabs(iti)) + enddo + if (surfarea) call sdarea(gamvec) +c if (lprn) then +c write (iout,*) "Matrix A and vector gamma" +c do i=1,dimen1 +c write (iout,'(i2,$)') i +c do j=1,dimen +c write (iout,'(f4.1,$)') A(i,j) +c enddo +c write (iout,'(f8.3)') gamvec(i) +c enddo +c endif + if (lprn) then + write (iout,*) "Vector gamvec" + do i=1,dimen1 + write (iout,'(i5,f10.5)') i, gamvec(i) + enddo + endif + +c The friction matrix + do k=1,dimen + do i=1,dimen + dtdi=0.0d0 + do j=1,dimen1 + dtdi=dtdi+A(j,k)*A(j,i)*gamvec(j) + enddo + fricmat(k,i)=dtdi + enddo + enddo + + if (lprn) then + write (iout,'(//a)') "Matrix fricmat" + call matout2(dimen,dimen,maxres2,maxres2,fricmat) + endif + if (lang.eq.2 .or. lang.eq.3) then +c Mass-scale the friction matrix if non-direct integration will be performed + do i=1,dimen + do j=1,dimen + Ginvfric(i,j)=0.0d0 + do k=1,dimen + do l=1,dimen + Ginvfric(i,j)=Ginvfric(i,j)+ + & Gsqrm(i,k)*Gsqrm(l,j)*fricmat(k,l) + enddo + enddo + enddo + enddo +c Diagonalize the friction matrix + ind=0 + do i=1,dimen + do j=1,i + ind=ind+1 + Ghalf(ind)=Ginvfric(i,j) + enddo + enddo + call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec, + & ierr,iwork) + if (lprn) then + write (iout,'(//2a)') "Eigenvectors and eigenvalues of the", + & " mass-scaled friction matrix" + call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam) + endif +c Precompute matrices for tinker stochastic integrator +#ifndef LANG0 + do i=1,dimen + do j=1,dimen + mt1(i,j)=0.0d0 + mt2(i,j)=0.0d0 + do k=1,dimen + mt1(i,j)=mt1(i,j)+fricvec(k,i)*gsqrm(k,j) + mt2(i,j)=mt2(i,j)+fricvec(k,i)*gsqrp(k,j) + enddo + mt3(j,i)=mt1(i,j) + enddo + enddo +#endif + else if (lang.eq.4) then +c Diagonalize the friction matrix + ind=0 + do i=1,dimen + do j=1,i + ind=ind+1 + Ghalf(ind)=fricmat(i,j) + enddo + enddo + call gldiag(maxres2,dimen,dimen,Ghalf,work,fricgam,fricvec, + & ierr,iwork) + if (lprn) then + write (iout,'(//2a)') "Eigenvectors and eigenvalues of the", + & " friction matrix" + call eigout(dimen,dimen,maxres2,maxres2,fricvec,fricgam) + endif +c Determine the number of zero eigenvalues of the friction matrix + nzero=max0(dimen-dimen1,0) +c do while (fricgam(nzero+1).le.1.0d-5 .and. nzero.lt.dimen) +c nzero=nzero+1 +c enddo + write (iout,*) "Number of zero eigenvalues:",nzero + do i=1,dimen + do j=1,dimen + fricmat(i,j)=0.0d0 + do k=nzero+1,dimen + fricmat(i,j)=fricmat(i,j) + & +fricvec(i,k)*fricvec(j,k)/fricgam(k) + enddo + enddo + enddo + if (lprn) then + write (iout,'(//a)') "Generalized inverse of fricmat" + call matout(dimen,dimen,maxres6,maxres6,fricmat) + endif + endif +#ifdef MPI + 10 continue + if (nfgtasks.gt.1) then + if (fg_rank.eq.0) then +c The matching BROADCAST for fg processors is called in ERGASTULUM +#ifdef MPI + time00=MPI_Wtime() +#else + time00=tcpu() +#endif + call MPI_Bcast(10,1,MPI_INTEGER,king,FG_COMM,IERROR) +#ifdef MPI + time_Bcast=time_Bcast+MPI_Wtime()-time00 +#else + time_Bcast=time_Bcast+tcpu()-time00 +#endif +c print *,"Processor",myrank, +c & " BROADCAST iorder in SETUP_FRICMAT" + endif +c licznik=licznik+1 +c write (iout,*) "setup_fricmat licznik",licznik +#ifdef MPI + time00=MPI_Wtime() +#else + time00=tcpu() +#endif +c 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) +#ifdef TIMING +#ifdef MPI + time_scatter=time_scatter+MPI_Wtime()-time00 + time_scatter_fmat=time_scatter_fmat+MPI_Wtime()-time00 +#else + time_scatter=time_scatter+tcpu()-time00 + time_scatter_fmat=time_scatter_fmat+tcpu()-time00 +#endif +#endif + do i=1,dimen + do j=1,2*my_ng_count + fricmat(j,i)=fcopy(i,j) + enddo + enddo +c write (iout,*) "My chunk of fricmat" +c call MATOUT2(my_ng_count,dimen,maxres2,maxres2,fcopy) + endif +#endif + return + end +c------------------------------------------------------------------------------- + subroutine sdarea(gamvec) +c +c Scale the friction coefficients according to solvent accessible surface areas +c Code adapted from TINKER +c AL 9/3/04 +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' +#ifndef LANG0 + include 'COMMON.LANGEVIN' +#else + include 'COMMON.LANGEVIN.lang0' +#endif + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision radius(maxres2),gamvec(maxres2) + parameter (twosix=1.122462048309372981d0) + logical lprn /.false./ +c +c determine new friction coefficients every few SD steps +c +c set the atomic radii to estimates of sigma values +c +c print *,"Entered sdarea" + probe = 0.0d0 + + do i=1,2*nres + radius(i)=0.0d0 + enddo +c Load peptide group radii + do i=nnt,nct-1 + radius(i)=pstok + enddo +c Load side chain radii + do i=nnt,nct + iti=itype(i) + radius(i+nres)=restok(iti) + enddo +c do i=1,2*nres +c write (iout,*) "i",i," radius",radius(i) +c enddo + do i = 1, 2*nres + radius(i) = radius(i) / twosix + if (radius(i) .ne. 0.0d0) radius(i) = radius(i) + probe + end do +c +c scale atomic friction coefficients by accessible area +c + if (lprn) write (iout,*) + & "Original gammas, surface areas, scaling factors, new gammas, ", + & "std's of stochastic forces" + ind=0 + do i=nnt,nct-1 + if (radius(i).gt.0.0d0) then + call surfatom (i,area,radius) + ratio = dmax1(area/(4.0d0*pi*radius(i)**2),1.0d-1) + if (lprn) write (iout,'(i5,3f10.5,$)') + & i,gamvec(ind+1),area,ratio + do j=1,3 + ind=ind+1 + gamvec(ind) = ratio * gamvec(ind) + enddo + stdforcp(i)=stdfp*dsqrt(gamvec(ind)) + if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcp(i) + endif + enddo + do i=nnt,nct + if (radius(i+nres).gt.0.0d0) then + call surfatom (i+nres,area,radius) + ratio = dmax1(area/(4.0d0*pi*radius(i+nres)**2),1.0d-1) + if (lprn) write (iout,'(i5,3f10.5,$)') + & i,gamvec(ind+1),area,ratio + do j=1,3 + ind=ind+1 + gamvec(ind) = ratio * gamvec(ind) + enddo + stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind)) + if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i) + endif + enddo + + return + end