chemshell_qchem / blob
1 subroutine nve_0
2 x (idnode,mxnode,natms,imcon,tstep,engke,weight,cell,
3 x xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz,buffer,stress)
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***********************************************************************
22 c
24 implicit real*8 (a-h,o-z)
25 include 'dl_params.inc'
27 dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
28 dimension vxx(mxatms),vyy(mxatms),vzz(mxatms)
29 dimension fxx(mxatms),fyy(mxatms),fzz(mxatms)
30 dimension weight(mxatms),buffer(mxbuff),cell(9)
31 dimension stress(9),stres1(9)
32 c
33 c set kinetic energy accumulator
35 engke=0.d0
36 #ifdef STRESS
37 c
38 c temporary stress tensor accumulators
39 do i = 1,9
40 stres1(i) = 0.d0
41 enddo
42 #endif
43 c
44 c block indices
46 iatm0 = (idnode*natms)/mxnode + 1
47 iatm1 = ((idnode+1)*natms)/mxnode
49 c
50 c move atoms by leapfrog algorithm
52 do i=iatm0,iatm1
53 c
54 c store old velocities
56 uxx=vxx(i)
57 uyy=vyy(i)
58 uzz=vzz(i)
59 c
60 c calculate new velocities
62 vxx(i)=vxx(i)+(tstep/weight(i))*fxx(i)
63 vyy(i)=vyy(i)+(tstep/weight(i))*fyy(i)
64 vzz(i)=vzz(i)+(tstep/weight(i))*fzz(i)
65 c
66 c update positions
68 xxx(i)=xxx(i)+tstep*vxx(i)
69 yyy(i)=yyy(i)+tstep*vyy(i)
70 zzz(i)=zzz(i)+tstep*vzz(i)
71 c
72 c current velocity
73 vxt = 0.5d0*(vxx(i)+uxx)
74 vyt = 0.5d0*(vyy(i)+uyy)
75 vzt = 0.5d0*(vzz(i)+uzz)
76 c
77 c calculate kinetic energy
78 engke=engke+0.5d0*weight(i)*(vxt*vxt+vyt*vyt+vzt*vzt)
80 #ifdef STRESS
81 c
82 c kinetic contribution to stress tensor
84 stres1(1) = stres1(1) + weight(i)*vxt*vxt
85 stres1(2) = stres1(2) + weight(i)*vxt*vyt
86 stres1(3) = stres1(3) + weight(i)*vxt*vzt
87 stres1(5) = stres1(5) + weight(i)*vyt*vyt
88 stres1(6) = stres1(6) + weight(i)*vyt*vzt
89 stres1(9) = stres1(9) + weight(i)*vzt*vzt
90 #endif
91 enddo
93 c
94 c periodic boundary condition
96 call images(imcon,idnode,mxnode,natms,cell,xxx,yyy,zzz)
98 c
99 c global exchange of configuration data
101 if(mxnode.gt.1)then
103 nbuff=mxbuff
104 call gdsum(engke,1,buffer)
105 call merge(idnode,mxnode,natms,nbuff,xxx,yyy,zzz,buffer)
106 call merge(idnode,mxnode,natms,nbuff,vxx,vyy,vzz,buffer)
108 endif
109 #ifdef STRESS
110 c
111 c complete stress tensor
113 if(mxnode.gt.1) call gdsum(stres1,9,buffer)
115 stress(1) = stress(1) + stres1(1)
116 stress(2) = stress(2) + stres1(2)
117 stress(3) = stress(3) + stres1(3)
118 stress(4) = stress(4) + stres1(2)
119 stress(5) = stress(5) + stres1(5)
120 stress(6) = stress(6) + stres1(6)
121 stress(7) = stress(7) + stres1(3)
122 stress(8) = stress(8) + stres1(6)
123 stress(9) = stress(9) + stres1(9)
124 #endif
125 return
126 end
