Skip to content

GitLab

Switch branch/tag
Find file Normal viewHistoryPermalink
kCalculateAmoebaCudaPmeDirectElectrostaticT2.h 7.94 KB
Newer Older
Mark Friedrichs's avatar
Optimzations and bug fix for mapping of torque to force for small molecules  
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 

static __device__ void SUB_METHOD_NAME( calculatePmeDirectElectrostaticPairIxnT2, _kernel)( 
                                        PmeDirectElectrostaticParticle& atomI, PmeDirectElectrostaticParticle& atomJ,
                                        const float4 delta, const float4 bn
#ifdef APPLY_SCALE
                                        , const float* scalingFactors
#endif
                                        ){

    float xr                    = delta.x;
    float yr                    = delta.y;
    float zr                    = delta.z;
    float rr1                   = delta.w;

    // set the permanent multipole and induced dipole values;

    float di1                   = atomI.labFrameDipole[0];
    float di2                   = atomI.labFrameDipole[1];
    float di3                   = atomI.labFrameDipole[2];

    float qi1                   = atomI.labFrameQuadrupole[0];
    float qi2                   = atomI.labFrameQuadrupole[1];
    float qi3                   = atomI.labFrameQuadrupole[2];
    float qi5                   = atomI.labFrameQuadrupole[3];
    float qi6                   = atomI.labFrameQuadrupole[4];
    //float qi9                   = atomI.labFrameQuadrupole[5];
    float qi9                   = -(atomI.labFrameQuadrupole[0] + atomI.labFrameQuadrupole[3]);

    float bn1                   = bn.x;
    float bn2                   = bn.y;
    float bn3                   = bn.z;

    // apply Thole polarization damping to scale factors

    float scale3                = 1.0f;
    float scale5                = 1.0f;
    float scale7                = 1.0f;

    float damp                  = atomI.damp*atomJ.damp;
    if( damp != 0.0f ){
        float pgamma  = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;
        float ratio   = 1.0f/(rr1*damp);
            damp      = -pgamma*ratio*ratio*ratio;
        if( damp > -50.0f ){
Mark Friedrichs's avatar
Renamed math functions in Cuda kernels to single-precision equivalent [cos() -> cosf()]  
Mark Friedrichs committed
45 
            float expdamp    = expf(damp);
Mark Friedrichs's avatar
Optimzations and bug fix for mapping of torque to force for small molecules  
Mark Friedrichs committed
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181 
            scale3           = 1.0f - expdamp;
            scale5           = 1.0f - (1.0f-damp)*expdamp;
            scale7           = 1.0f - (1.0f-damp+0.6f*damp*damp)*expdamp;
        }
    }


    float rr3                   = rr1*rr1*rr1;
#ifdef APPLY_SCALE
    float dsc3                  = rr3*(1.0f - scale3*scalingFactors[DScaleIndex]);
    float dsc5                  = (3.0f*rr3*rr1*rr1)* (1.0f - scale5*scalingFactors[DScaleIndex]);
    float dsc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7*scalingFactors[DScaleIndex]);

    float psc3                  = rr3*(1.0f - scale3*scalingFactors[PScaleIndex]);
    float psc5                  = (3.0f*rr3*rr1*rr1)*(1.0f - scale5*scalingFactors[PScaleIndex]);
    float psc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7*scalingFactors[PScaleIndex]);
#else
    float psc3                  = rr3*(1.0f - scale3);
    float psc5                  = (3.0f*rr3*rr1*rr1)*(1.0f - scale5);
    float psc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7);
#endif

    float prefactor1            = 0.5f*(psc3 - bn1);
#ifdef APPLY_SCALE
    float prefactor2            = 0.5f*(dsc3 - bn1);
#endif

    float dixuk1                = di2*atomJ.inducedDipole[2]  - di3*atomJ.inducedDipole[1];
    float dixukp1               = di2*atomJ.inducedDipoleP[2] - di3*atomJ.inducedDipoleP[1];

#ifdef APPLY_SCALE
    float ttm2i1                = prefactor1*dixuk1 + prefactor2*dixukp1;
#else
    float ttm2i1                = prefactor1*(dixuk1 + dixukp1);
#endif

    float dixuk2                = di3*atomJ.inducedDipole[0]  - di1*atomJ.inducedDipole[2];
    float dixukp2               = di3*atomJ.inducedDipoleP[0] - di1*atomJ.inducedDipoleP[2];

#ifdef APPLY_SCALE
    float ttm2i2                = prefactor1*dixuk2 + prefactor2*dixukp2;
#else
    float ttm2i2                = prefactor1*(dixuk2 + dixukp2);
#endif

    float dixuk3                = di1*atomJ.inducedDipole[1]  - di2*atomJ.inducedDipole[0];
    float dixukp3               = di1*atomJ.inducedDipoleP[1] - di2*atomJ.inducedDipoleP[0];
#ifdef APPLY_SCALE
    float ttm2i3                = prefactor1*dixuk3 + prefactor2*dixukp3;
#else
    float ttm2i3                = prefactor1*(dixuk3 + dixukp3);
#endif

    float sci4                  = atomJ.inducedDipole[0]*xr  + atomJ.inducedDipole[1]*yr  + atomJ.inducedDipole[2]*zr;
    float scip4                 = atomJ.inducedDipoleP[0]*xr + atomJ.inducedDipoleP[1]*yr + atomJ.inducedDipoleP[2]*zr;
    float gti2                  = bn2*(sci4+scip4);
#ifdef APPLY_SCALE
    float gtri2                 = (sci4*psc5+scip4*dsc5);
#else
    float gtri2                 = psc5*(sci4+scip4);
#endif
    prefactor1                  = 0.5f*(gti2 - gtri2);

    ttm2i1                     += prefactor1*( di2*zr - di3*yr );
    ttm2i2                     += prefactor1*( di3*xr - di1*zr );
    ttm2i3                     += prefactor1*( di1*yr - di2*xr );

    float qir1                  = qi1*xr + qi2*yr + qi3*zr;
    float qir2                  = qi2*xr + qi5*yr + qi6*zr;
    float qir3                  = qi3*xr + qi6*yr + qi9*zr;

#ifdef APPLY_SCALE
    prefactor1                  = sci4*psc7 + scip4*dsc7 - bn3*(sci4+scip4);
#else
    prefactor1                  = psc7*(sci4+scip4) - bn3*(sci4+scip4);
#endif
    ttm2i1                     += prefactor1*( yr*qir3 - zr*qir2 );
    ttm2i2                     += prefactor1*( zr*qir1 - xr*qir3 );
    ttm2i3                     += prefactor1*( xr*qir2 - yr*qir1 );

    float qiuk1                 = qi1*atomJ.inducedDipole[0]  + qi2*atomJ.inducedDipole[1]  + qi3*atomJ.inducedDipole[2];
    float qiuk2                 = qi2*atomJ.inducedDipole[0]  + qi5*atomJ.inducedDipole[1]  + qi6*atomJ.inducedDipole[2];
    float qiuk3                 = qi3*atomJ.inducedDipole[0]  + qi6*atomJ.inducedDipole[1]  + qi9*atomJ.inducedDipole[2];

    float qiukp1                = qi1*atomJ.inducedDipoleP[0] + qi2*atomJ.inducedDipoleP[1] + qi3*atomJ.inducedDipoleP[2];
    float qiukp2                = qi2*atomJ.inducedDipoleP[0] + qi5*atomJ.inducedDipoleP[1] + qi6*atomJ.inducedDipoleP[2];
    float qiukp3                = qi3*atomJ.inducedDipoleP[0] + qi6*atomJ.inducedDipoleP[1] + qi9*atomJ.inducedDipoleP[2];

    prefactor1                  = (bn2 - psc5);
#ifdef APPLY_SCALE
    prefactor2                  = (bn2 - dsc5);
#endif
    float ukxqir1               = atomJ.inducedDipole[1]*qir3  - atomJ.inducedDipole[2]*qir2;
    float ukxqirp1              = atomJ.inducedDipoleP[1]*qir3 - atomJ.inducedDipoleP[2]*qir2;
    float rxqiuk1               = yr*qiuk3  - zr*qiuk2;
    float rxqiukp1              = yr*qiukp3 - zr*qiukp2;

#ifdef APPLY_SCALE
    ttm2i1                     += prefactor1*(ukxqir1 + rxqiuk1) + prefactor2*(ukxqirp1 + rxqiukp1);
#else
    ttm2i1                     += prefactor1*( ukxqir1 + rxqiuk1 + ukxqirp1 + rxqiukp1 );
#endif

    float ukxqir2               = atomJ.inducedDipole[2]*qir1  - atomJ.inducedDipole[0]*qir3;
    float ukxqirp2              = atomJ.inducedDipoleP[2]*qir1 - atomJ.inducedDipoleP[0]*qir3;
    float rxqiuk2               = zr*qiuk1  - xr*qiuk3;
    float rxqiukp2              = zr*qiukp1 - xr*qiukp3;
#ifdef APPLY_SCALE
    ttm2i2                     += prefactor1*(ukxqir2 + rxqiuk2) + prefactor2*(ukxqirp2 + rxqiukp2);
#else
    ttm2i2                     += prefactor1*( ukxqir2 + rxqiuk2 + ukxqirp2 + rxqiukp2 );
#endif

    float ukxqir3               = atomJ.inducedDipole[0]*qir2  - atomJ.inducedDipole[1]*qir1;
    float ukxqirp3              = atomJ.inducedDipoleP[0]*qir2 - atomJ.inducedDipoleP[1]*qir1;
    float rxqiuk3               = xr*qiuk2  - yr*qiuk1;
    float rxqiukp3              = xr*qiukp2 - yr*qiukp1;
#ifdef APPLY_SCALE
    ttm2i3                     += prefactor1*(ukxqir3 + rxqiuk3) + prefactor2*(ukxqirp3 + rxqiukp3);
#else
    ttm2i3                     += prefactor1*(ukxqir3 + rxqiuk3 + ukxqirp3 + rxqiukp3 );
#endif

    atomI.torque[0]            += ttm2i1;
    atomI.torque[1]            += ttm2i2;
    atomI.torque[2]            += ttm2i3;

/*
    torque[0]            += ttm2i1;
    torque[1]            += ttm2i2;
    torque[2]            += ttm2i3;
*/

    return;

}