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.21) 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.21) 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.21) 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.21) 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.21) 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 do i=nnt,nct ind=ind+1 ii = ind+m iti=itype(i) gamvec(ii)=gamsc(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_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