mpi2/forces.cpp at master · jonqdoe/mpi2 · GitHub

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
#include "globals.h"
void charge_grid( void ) ;
void bonds( void ) ;
void charge_forces( void ) ;
void communicate_ghost_forces( void ) ;
void angles( void ) ;
void particle_orientations( void ) ;
void calc_S_conv_uG(void);
void ms_forces( void ) ;


void forces() {

  int i,j, m, gind, t1, t2, id ;


  charge_grid() ;

  ///////////////////////////////////////////////
  // Reset the particle forces and grid grad w //
  ///////////////////////////////////////////////
  for ( i=0 ; i<ns_loc ; i++ ) {
    id = my_inds[i] ;
    for ( j=0 ; j<Dim ; j++ )
      f[id][j] = 0.0 ;
  }
  for ( i=0 ; i<total_ghost ; i++ ) {
    id = ghost_inds[i] ;
    for ( j=0 ; j<Dim ; j++ )
      f[id][j] = 0.0 ;
  }


  for ( j=0 ; j<ntypes ; j++ )
    Components[j].ZeroGradient() ;


  /////////////////////////////////////
  // Calculate the grid-based forces //
  /////////////////////////////////////
  for ( j=0 ; j<n_gaussian_pairstyles ; j++ )
    Gauss[j].CalcAll() ;

  for ( j=0 ; j<n_gausserfc_pairstyles ; j++ )
    GaussErfc[j].CalcAll() ;

  for ( j=0 ; j<n_erfc2_pairstyles ; j++ )
    Erfc2[j].CalcAll() ;


  //////////////////////////////////////////////////////
  // Accumulate the nonbonded forces on each particle //
  //////////////////////////////////////////////////////
  for ( i=0 ; i<ns_loc+total_ghost ; i++ ) {

    if ( i < ns_loc )
      id = my_inds[i] ;
    else
      id = ghost_inds[ i-ns_loc ] ;


    for ( m=0 ; m < grid_per_partic ; m++ ) {
      // Note that in charge_grid(), grid_inds is populated with
      // indices in the range [0,ML].
      // if gind == -1, then grid location is skipped.
      gind = grid_inds[ id ][ m ] ;

      if ( gind == -1 )
        continue ;

      for ( j=0 ; j<Dim ; j++ ) {
        // This division by rho(r) might be better off pre-calculated
        f[id][j] += Components[ tp[id] ].force[j][gind] / Components[ tp[id] ].rho[gind] ;
      }
    }

    for ( j=0 ; j<Dim ; j++ )
      f[id][j] *= gvol ;
  }


  ////////////////////////////
  // Call the bonded forces //
  ////////////////////////////
  if ( n_total_bonds > 0 )
    bonds() ;

  if ( n_total_angles > 0 )
    angles() ;


  if ( nprocs > 1 )
    communicate_ghost_forces() ;

}