Unverified Commit f380a03a authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1237 from akohlmey/compute-mop-bugfix 

Make compute stress/mop and stress/mop/profile compatible with per-atom masses
parents 1f210a24 7b68655c

src/USER-MISC/compute_stress_mop.cpp

+24 −11
Original line number Diff line number Diff line
@@ -252,6 +252,7 @@ void ComputeStressMop::compute_pairs() 
  int *ilist,*jlist,*numneigh,**firstneigh;



  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;
@@ -388,9 +389,15 @@ void ComputeStressMop::compute_pairs() 
          fi[2] = atom->f[i][2];



          //coordinates at t-dt (based on Velocity-Verlet alg.)

          if (rmass) {

            xj[0] = xi[0]-vi[0]*dt+fi[0]/2/rmass[i]*dt*dt*ftm2v;

            xj[1] = xi[1]-vi[1]*dt+fi[1]/2/rmass[i]*dt*dt*ftm2v;

            xj[2] = xi[2]-vi[2]*dt+fi[2]/2/rmass[i]*dt*dt*ftm2v;

          } else {

            xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;

            xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;

            xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;

          }



          // because LAMMPS does not put atoms back in the box

          // at each timestep, must check atoms going through the
@@ -406,9 +413,15 @@ void ComputeStressMop::compute_pairs() 


            //approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)

            double vcross[3];

            if (rmass) {

              vcross[0] = vi[0]-fi[0]/rmass[i]/2*ftm2v*dt;

              vcross[1] = vi[1]-fi[1]/rmass[i]/2*ftm2v*dt;

              vcross[2] = vi[2]-fi[2]/rmass[i]/2*ftm2v*dt;

            } else {

              vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;

              vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;

              vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;

            }



            values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;

            values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;

src/USER-MISC/compute_stress_mop_profile.cpp

+35 −17
Original line number Diff line number Diff line
@@ -260,6 +260,7 @@ void ComputeStressMopProfile::compute_pairs() 
  double pos,pos1;



  double *mass = atom->mass;

  double *rmass = atom->rmass;

  int *type = atom->type;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;
@@ -411,9 +412,15 @@ void ComputeStressMopProfile::compute_pairs() 
          fi[2] = atom->f[i][2];



          //coordinates at t-dt (based on Velocity-Verlet alg.)

          if (rmass) {

            xj[0] = xi[0]-vi[0]*dt+fi[0]/2/rmass[i]*dt*dt*ftm2v;

            xj[1] = xi[1]-vi[1]*dt+fi[1]/2/rmass[i]*dt*dt*ftm2v;

            xj[2] = xi[2]-vi[2]*dt+fi[2]/2/rmass[i]*dt*dt*ftm2v;

          } else {

            xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v;

            xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v;

            xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v;

          }



          for (ibin=0;ibin<nbins;ibin++) {

            pos = coord[ibin][0];
@@ -424,6 +431,16 @@ void ComputeStressMopProfile::compute_pairs() 
              sgn = copysign(1.0,vi[dir]);



              //approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.)

              if (rmass) {

                double vcross[3];

                vcross[0] = vi[0]-fi[0]/rmass[i]/2*ftm2v*dt;

                vcross[1] = vi[1]-fi[1]/rmass[i]/2*ftm2v*dt;

                vcross[2] = vi[2]-fi[2]/rmass[i]/2*ftm2v*dt;



                values_local[ibin][m] += rmass[i]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;

                values_local[ibin][m+1] += rmass[i]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;

                values_local[ibin][m+2] += rmass[i]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;

              } else {

                double vcross[3];

                vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;

                vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
@@ -437,6 +454,7 @@ void ComputeStressMopProfile::compute_pairs() 
          }

        }

      }

    }

    m+=3;

  }

}

src/compute_reduce_chunk.h

+1 −1

File changed.

Contains only whitespace changes.

src/compute_adf.h

+1 −1

File changed.

Contains only whitespace changes.

CRA Git | Maintained and supported by SUSTech CRA and CCSE