chemshell_qchem / blob
1 subroutine friction_0
2 x (idnode,mxnode,natms,imcon,tstep,engke,weight,cell,
3 x xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz,buffer,stress,fric)
5 c
6 c***********************************************************************
7 c
8 c dl_poly subroutine for integrating newtonian equations of
9 c motion in molecular dynamics - verlet leapfrog.
10 c
11 c parallel replicated data version
12 c
13 c copyright - daresbury laboratory 1992
14 c author - w. smith march 1992.
15 c
16 c $Author: thiels $
17 c $Date: 2004-08-19 13:51:33 +0200 (Thu, 19. Aug 2004) $
18 c $Revision: 2490 $
19 c $State$
20 c
21 c FRICTION DYNAMICS
22 c Propagate with a constant friction proportional to the velocity.
23 c The friction parameter fric is defined as
24 c fric = friction coefficient * Dt / 2
25 c
26 c Hans Martin Senn (hms), MPI fuer Kohlenforschung, Nov. 2002
27 c***********************************************************************
28 c
30 implicit real*8 (a-h,o-z)
31 include 'dl_params.inc'
33 dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
34 dimension vxx(mxatms),vyy(mxatms),vzz(mxatms)
35 dimension fxx(mxatms),fyy(mxatms),fzz(mxatms)
36 dimension weight(mxatms),buffer(mxbuff),cell(9)
37 dimension stress(9),stres1(9)
38 c
39 c set kinetic energy accumulator
41 engke=0.d0
42 #ifdef STRESS
43 c
44 c temporary stress tensor accumulators
45 do i = 1,9
46 stres1(i) = 0.d0
47 enddo
48 #endif
49 c
50 c block indices
52 iatm0 = (idnode*natms)/mxnode + 1
53 iatm1 = ((idnode+1)*natms)/mxnode
55 c
56 c move atoms by leapfrog algorithm
57 svar = 1.d0/(1.d0 + fric)
59 do i=iatm0,iatm1
60 c
61 c store old velocities
63 uxx=vxx(i)
64 uyy=vyy(i)
65 uzz=vzz(i)
66 c
67 c calculate new velocities
68 svar2 = tstep/weight(i)
69 vxx(i)=svar * ((1.d0 - fric)*vxx(i) + svar2*fxx(i))
70 vyy(i)=svar * ((1.d0 - fric)*vyy(i) + svar2*fyy(i))
71 vzz(i)=svar * ((1.d0 - fric)*vzz(i) + svar2*fzz(i))
72 c
73 c update positions
75 xxx(i)=xxx(i)+tstep*vxx(i)
76 yyy(i)=yyy(i)+tstep*vyy(i)
77 zzz(i)=zzz(i)+tstep*vzz(i)
78 c
79 c current velocity
80 vxt = 0.5d0*(vxx(i)+uxx)
81 vyt = 0.5d0*(vyy(i)+uyy)
82 vzt = 0.5d0*(vzz(i)+uzz)
83 c
84 c calculate kinetic energy
85 engke=engke+0.5d0*weight(i)*(vxt*vxt+vyt*vyt+vzt*vzt)
87 #ifdef STRESS
88 c
89 c kinetic contribution to stress tensor
91 stres1(1) = stres1(1) + weight(i)*vxt*vxt
92 stres1(2) = stres1(2) + weight(i)*vxt*vyt
93 stres1(3) = stres1(3) + weight(i)*vxt*vzt
94 stres1(5) = stres1(5) + weight(i)*vyt*vyt
95 stres1(6) = stres1(6) + weight(i)*vyt*vzt
96 stres1(9) = stres1(9) + weight(i)*vzt*vzt
97 #endif
98 enddo
100 c
101 c periodic boundary condition
103 call images(imcon,idnode,mxnode,natms,cell,xxx,yyy,zzz)
105 c
106 c global exchange of configuration data
108 if(mxnode.gt.1)then
110 nbuff=mxbuff
111 call gdsum(engke,1,buffer)
112 call merge(idnode,mxnode,natms,nbuff,xxx,yyy,zzz,buffer)
113 call merge(idnode,mxnode,natms,nbuff,vxx,vyy,vzz,buffer)
115 endif
116 #ifdef STRESS
117 c
118 c complete stress tensor
120 if(mxnode.gt.1) call gdsum(stres1,9,buffer)
122 stress(1) = stress(1) + stres1(1)
123 stress(2) = stress(2) + stres1(2)
124 stress(3) = stress(3) + stres1(3)
125 stress(4) = stress(4) + stres1(2)
126 stress(5) = stress(5) + stres1(5)
127 stress(6) = stress(6) + stres1(6)
128 stress(7) = stress(7) + stres1(3)
129 stress(8) = stress(8) + stres1(6)
130 stress(9) = stress(9) + stres1(9)
131 #endif
132 return
133 end
