subroutine friction_force implicit none include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #endif include 'COMMON.IOUNITS' #ifdef FIVEDIAG integer iposc,ichain,n,innt,inct double precision v_work(3,maxres2),vvec(maxres2),rs(maxres2) #else double precision gamvec(MAXRES6) common /syfek/ gamvec double precision vv(3),vvtot(3,maxres),v_work(MAXRES6), & ginvfric(maxres2,maxres2) common /przechowalnia/ ginvfric #endif integer i,j,k,ind logical lprn /.false./, checkmode /.false./ #ifdef FIVEDIAG c Here accelerations due to friction forces are computed right after forces. d_t_work=0.0d0 do j=1,3 v_work(j,1)=d_t(j,0) v_work(j,nnt)=d_t(j,0) enddo do i=nnt+1,nct do j=1,3 v_work(j,i)=v_work(j,i-1)+d_t(j,i-1) enddo enddo do i=nnt,nct if (iabs(itype(i)).ne.10 .and. iabs(itype(i)).ne.ntyp1) then do j=1,3 v_work(j,i+nres)=v_work(j,i)+d_t(j,i+nres) enddo endif enddo #ifdef DEBUG write (iout,*) "v_work" do i=1,2*nres write (iout,'(i5,3f10.5)') i,(v_work(j,i),j=1,3) enddo #endif do j=1,3 ind=0 do ichain=1,nchain n=dimen_chain(ichain) iposc=iposd_chain(ichain) c write (iout,*) "friction_force j",j," ichain",ichain, c & " n",n," iposc",iposc,iposc+n-1 innt=chain_border(1,ichain) inct=chain_border(2,ichain) c diagnostics c innt=chain_border(1,1) c inct=chain_border(2,1) do i=innt,inct vvec(ind+1)=v_work(j,i) ind=ind+1 if (iabs(itype(i)).ne.10) then vvec(ind+1)=v_work(j,i+nres) ind=ind+1 endif enddo #ifdef DEBUG write (iout,*) "vvec ind",ind," n",n write (iout,'(f10.5)') (vvec(i),i=iposc,ind) #endif c write (iout,*) "chain",i," ind",ind," n",n call fivediagmult(n,DMfric(iposc),DU1fric(iposc), & DU2fric(iposc),vvec(iposc),rs) #ifdef DEBUG write (iout,*) "rs" write (iout,'(f10.5)') (rs(i),i=1,n) #endif do i=iposc,iposc+n-1 c write (iout,*) "ichain",ichain," i",i," j",j, c & "index",3*(i-1)+j,"rs",rs(i-iposc+1) fric_work(3*(i-1)+j)=-rs(i-iposc+1) enddo enddo enddo #ifdef DEBUG write (iout,*) "Vector fric_work dimen3",dimen3 write (iout,'(3f10.5)') (fric_work(j),j=1,dimen3) #endif #else 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 #endif return end c----------------------------------------------------- subroutine stochastic_force(stochforcvec) implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' double precision time00 #endif include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif include 'COMMON.TIME1' #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #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./ integer i,j,ind,ii,iti double precision anorm_distr #ifdef FIVEDIAG integer ichain,innt,inct,iposc #endif 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 #ifdef FIVEDIAG if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle #endif 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 DEBUG write (iout,*) "Stochastic forces on sites" do i=1,nres write (iout,'(i5,2(3f10.5,5x))') i,(force(j,i),j=1,3), & (force(j,i+nres),j=1,3) enddo #endif #ifdef MPI time_fsample=time_fsample+MPI_Wtime()-time00 #else time_fsample=time_fsample+tcpu()-time00 #endif #ifdef FIVEDIAG ind=0 do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) iposc=iposd_chain(ichain) c for debugging only c innt=chain_border(1,1) c inct=chain_border(2,1) c iposc=iposd_chain(1) c write (iout,*)"stochastic_force ichain=",ichain," innt",innt, c & " inct",inct," iposc",iposc do j=1,3 stochforcvec(ind+j)=0.5d0*force(j,innt) enddo if (iabs(itype(innt)).eq.10) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,innt+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,innt+nres) enddo ind=ind+3 endif do i=innt+1,inct-1 do j=1,3 stochforcvec(ind+j)=0.5d0*(force(j,i)+force(j,i-1)) enddo if (iabs(itype(i)).eq.10) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,i+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,i+nres) enddo ind=ind+3 endif enddo do j=1,3 stochforcvec(ind+j)=0.5d0*force(j,inct-1) enddo if (iabs(itype(inct)).eq.10) then do j=1,3 stochforcvec(ind+j)=stochforcvec(ind+j)+force(j,inct+nres) enddo ind=ind+3 else ind=ind+3 do j=1,3 stochforcvec(ind+j)=force(j,inct+nres) enddo ind=ind+3 endif c write (iout,*) "chain",ichain," ind",ind enddo #ifdef DEBUG write (iout,*) "stochforcvec" write (iout,'(3f10.5)') (stochforcvec(j),j=1,ind) #endif #else 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 #endif return end c------------------------------------------------------------------ subroutine setup_fricmat implicit none #ifdef MPI include 'mpif.h' integer ierr double precision time00 #endif include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' #ifdef FIVEDIAG include 'COMMON.LAGRANGE.5diag' #else include 'COMMON.LAGRANGE' #endif include 'COMMON.SETUP' include 'COMMON.TIME1' c integer licznik /0/ c save licznik #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #endif include 'COMMON.IOUNITS' integer IERROR integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct integer ichain,nind logical lprn /.false./ double precision dtdi,gamvec(MAXRES2) common /syfek/ gamvec #ifndef FIVEDIAG double precision ginvfric(maxres2,maxres2),Ghalf(mmaxres2), & fcopy(maxres2,maxres2) double precision work(8*maxres2) integer iwork(maxres2) common /przechowalnia/ ginvfric,Ghalf,fcopy #endif #ifdef MPI if (fg_rank.ne.king) goto 10 #endif 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 if (lprn) write (iout,*) "m",m ind=0 C 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) if (lprn) then write (iout,*) "Vector gamvec ii",ii do i=1,ii write (iout,'(i5,f10.5)') i, gamvec(i) enddo endif #ifdef FIVEDIAG DMfric=0.0d0 DU1fric=0.0d0 DU2fric=0.0d0 ind=1 do ichain=1,nchain innt=chain_border(1,ichain) inct=chain_border(2,ichain) c write (iout,*) "ichain",ichain," innt",innt," inct",inct c DMfric part DMfric(ind)=gamvec(innt-nnt+1)/4 if (iabs(itype(innt)).eq.10) then DMfric(ind)=DMfric(ind)+gamvec(m+innt-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(m+innt-nnt+1) ind=ind+2 endif c write (iout,*) "DMfric init ind",ind c DMfric do i=innt+1,inct-1 DMfric(ind)=gamvec(i-nnt+1)/2 if (iabs(itype(i)).eq.10) then DMfric(ind)=DMfric(ind)+gamvec(m+i-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(m+i-nnt+1) ind=ind+2 endif enddo c write (iout,*) "DMfric endloop ind",ind if (inct.gt.innt) then DMfric(ind)=gamvec(inct-1-nnt+1)/4 if (iabs(itype(inct)).eq.10) then DMfric(ind)=DMfric(ind)+gamvec(inct+m-nnt+1) ind=ind+1 else DMfric(ind+1)=gamvec(inct+m-nnt+1) ind=ind+2 endif endif c write (iout,*) "DMfric end ind",ind enddo c DU1fric part do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct if (iabs(itype(i)).ne.10) then ind=ind+2 else DU1fric(ind)=gamvec(i-nnt+1)/4 ind=ind+1 endif enddo enddo c DU2fric part do ichain=1,nchain ind=iposd_chain(ichain) innt=chain_border(1,ichain) inct=chain_border(2,ichain) do i=innt,inct-1 if (iabs(itype(i)).ne.10) then DU2fric(ind)=gamvec(i-nnt+1)/4 DU2fric(ind+1)=0.0d0 ind=ind+2 else DU2fric(ind)=0.0d0 ind=ind+1 endif enddo enddo if (lprn) then write(iout,*)"The upper part of the five-diagonal friction matrix" do ichain=1,nchain write (iout,'(a,i5)') 'Chain',ichain innt=iposd_chain(ichain) inct=iposd_chain(ichain)+dimen_chain(ichain)-1 do i=innt,inct if (i.lt.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i), & DU2fric(i) else if (i.eq.inct-1) then write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i),DU1fric(i) else write (iout,'(2i3,3f10.5)') i,i-innt+1,DMfric(i) endif enddo enddo endif 10 continue #else 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 #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 none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.VAR' include 'COMMON.MD' #ifndef LANG0 include 'COMMON.LANGEVIN' #else #ifdef FIVEDIAG include 'COMMON.LANGEVIN.lang0.5diag' #else include 'COMMON.LANGEVIN.lang0' #endif #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) double precision twosix parameter (twosix=1.122462048309372981d0) logical lprn /.false./ integer i,j,iti,ind double precision probe,area,ratio 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) if (iti.ne.ntyp1) 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 if (itype(i).ne.ntyp1) & stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamvec(ind)) if (lprn) write (iout,'(2f10.5)') gamvec(ind),stdforcsc(i) endif enddo return end