kupdatemd.br 3.03 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

/****************************************************************
* This file is part of the gpu acceleration library for gromacs.
* Author: V. Vishal
* Copyright (C) Pande Group, Stanford, 2006
*****************************************************************/

kernel void kupdate_md1_berendsen( 
		float dt,
		float3 lg, //Berendsen coupling, assuming only one group
		float3 uold, //Mean velocity from previous step
		float4 posq<>, 
		float3 v<>, 
		float3 f<>,
		float invmass<>,
		out float3 vnew<>, 
		out float4 posqp<> )
{
	float3 vb;
	float3 one;
	posqp = posq;

	one = float3( 1.0f, 1.0f, 1.0f );

	vb = ( one - lg ) * uold;
	
	vnew = lg* ( v + f * invmass * dt ) + vb;
 	
	posqp.xyz += vnew * dt;
}

//Nose-Hoover / Parinello Rahman
kernel void kupdate_md1_extended (
		float dt,
		float3 lg, 
		float xi,
		float3 M0, 
		float3 M1,
		float3 M2,
		float3 uold,
		float4 posq<>, 
		float3 v<>, 
		float3 f<>,
		float invmass<>,
		out float3 vnew<>, 
		out float4 posqp<> )
{
	float3 vrel, vnrel;
	float3 vtrans;

	vrel  = v - uold;
	vtrans = float3( dot(M0, vrel), dot(M1, vrel), dot(M2, vrel) );
	vrel += dt * ( invmass * f - xi * vrel - vtrans );
	vnew = uold + lg * vrel;
	posqp = posq;
	posqp.xyz += vnew * dt;
}

		
kernel void kupdate_md2(
		float dtinv, //1/dt
		float4 posqp<>, //positions after constraints
		float4 posq<>,  //positions before update
		out float3 vnew<>, //Corrected velocities
		out float4 posqnew<> //equal to posqp, avoids an extra call to copy
		)
{
	posqnew = posqp;
	vnew = ( posqp - posq ) * dtinv;
}

/* Version that handles terms of order dt more carefully
 * Update1 computes deltax's rather than x + deltax 
 * Corresponding shake implementation modifies the delta's
 * Update2 adds the constrained deltas to x.
 * */

kernel void kupdate_md1_berendsen_fix1( 
		float dt,
		float3 lg, //Berendsen coupling, assuming only one group
		float3 uold, //Mean velocity from previous step
		float4 posq<>, 
		float3 v<>, 
		float3 f<>,
		float invmass<>,
		out float3 vnew<>, 
		out float4 posqp<> )
{
	float3 vb;
	float3 one;

	one = float3( 1.0f, 1.0f, 1.0f );

	vb = ( one - lg ) * uold;
	
	vnew = lg* ( v + f * invmass * dt ) + vb;
 	
	posqp.xyz = vnew * dt;
	posqp.w = posq.w;
}

//Nose-Hoover / Parinello Rahman
kernel void kupdate_md1_extended_fix1 (
		float dt,
		float3 lg, 
		float xi,
		float3 M0, 
		float3 M1,
		float3 M2,
		float3 uold,
		float4 posq<>, 
		float3 v<>, 
		float3 f<>,
		float invmass<>,
		out float3 vnew<>, 
		out float4 posqp<> )
{
	float3 vrel, vnrel;
	float3 vtrans;

	vrel  = v - uold;
	vtrans = float3( dot(M0, vrel), dot(M1, vrel), dot(M2, vrel) );
	vrel += dt * ( invmass * f - xi * vrel - vtrans );
	vnew = uold + lg * vrel;
	posqp.w = posq.w;
	posqp.xyz = vnew * dt;
}

kernel void kupdate_md2_fix1(
		float dtinv, //1/dt
		float4 posqp<>, //positions after constraints
		float4 posq<>,  //positions before update
		out float3 vnew<>, //Corrected velocities
		out float4 posqnew<> //equal to posqp, avoids an extra call to copy
		)
{
	posqnew.xyz = posq.xyz + posqp.xyz;
	posqnew.w   = posq.w;
	vnew = posqp.xyz * dtinv;
}