(bc155995) · Commits · 郑智淋 / lammps

src/compute_angle_local.cpp

+31 −73

Original line number	Diff line number	Diff line
		@@ -25,8 +25,6 @@

		using namespace LAMMPS_NS;

		enum{THETA,ENERGY};

		#define DELTA 10000

		/* ---------------------------------------------------------------------- */
		@@ -44,28 +42,19 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS lmp, int narg, char *arg) :
		if (nvalues == 1) size_local_cols = 0;
		else size_local_cols = nvalues;

		which = new int[nvalues];
		pack_choice = new FnPtrPack[nvalues];

		tflag = eflag = 0;
		tflag = eflag = -1;
		nvalues = 0;

		int i;
		for (int iarg = 3; iarg < narg; iarg++) {
		i = iarg-3;

		if (strcmp(arg[iarg],"theta") == 0) {
		tflag = 1;
		which[i] = THETA;
		pack_choice[i] = &ComputeAngleLocal::pack_theta;
		} else if (strcmp(arg[iarg],"energy") == 0) {
		eflag = 1;
		which[i] = ENERGY;
		pack_choice[i] = &ComputeAngleLocal::pack_energy;
		} else error->all("Invalid keyword in compute angle/local command");
		if (strcmp(arg[iarg],"theta") == 0) tflag = nvalues++;
		else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
		else error->all("Invalid keyword in compute angle/local command");
		}

		nmax = 0;
		theta = energy = NULL;
		vector = NULL;
		array = NULL;
		}

		@@ -73,10 +62,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS lmp, int narg, char *arg) :

		ComputeAngleLocal::~ComputeAngleLocal()
		{
		delete [] which;
		delete [] pack_choice;
		memory->sfree(theta);
		memory->sfree(energy);
		memory->sfree(vector);
		memory->destroy_2d_double_array(array);
		}

		@@ -106,14 +92,6 @@ void ComputeAngleLocal::compute_local()
		if (ncount > nmax) reallocate(ncount);
		size_local_rows = ncount;
		ncount = compute_angles(1);

		// fill array with theta/energy values

		if (nvalues > 1) {
		if (array) buf = array[0];
		for (int n = 0; n < nvalues; n++)
		(this->*pack_choice[n])(n);
		}
		}

		/* ----------------------------------------------------------------------
		@@ -129,9 +107,10 @@ void ComputeAngleLocal::compute_local()

		int ComputeAngleLocal::compute_angles(int flag)
		{
		int i,atom1,atom2,atom3;
		int i,m,n,atom1,atom2,atom3;
		double delx1,dely1,delz1,delx2,dely2,delz2;
		double rsq1,rsq2,r1,r2,c;
		double tbuf,ebuf;

		double **x = atom->x;
		int *num_angle = atom->num_angle;
		@@ -143,10 +122,20 @@ int ComputeAngleLocal::compute_angles(int flag)
		int *mask = atom->mask;
		int nlocal = atom->nlocal;

		if (flag) {
		if (nvalues == 1) {
		if (tflag >= 0) tbuf = vector;
		if (eflag >= 0) ebuf = vector;
		} else {
		if (tflag >= 0) tbuf = &array[0][tflag];
		if (eflag >= 0) ebuf = &array[0][eflag];
		}
		}

		Angle *angle = force->angle;
		double PI = 4.0*atan(1.0);

		int m = 0;
		m = n = 0;
		for (atom2 = 0; atom2 < nlocal; atom2++) {
		if (!(mask[atom2] & groupbit)) continue;
		for (i = 0; i < num_angle[atom2]; i++) {
		@@ -158,7 +147,7 @@ int ComputeAngleLocal::compute_angles(int flag)
		if (angle_type[atom2][i] == 0) continue;

		if (flag) {
		if (tflag) {
		if (tflag >= 0) {
		delx1 = x[atom1][0] - x[atom2][0];
		dely1 = x[atom1][1] - x[atom2][1];
		delz1 = x[atom1][2] - x[atom2][2];
		@@ -181,14 +170,15 @@ int ComputeAngleLocal::compute_angles(int flag)
		c /= r1*r2;
		if (c > 1.0) c = 1.0;
		if (c < -1.0) c = -1.0;
		theta[m] = 180.0*acos(c)/PI;
		tbuf[n] = 180.0*acos(c)/PI;
		}

		if (eflag) {
		if (eflag >= 0) {
		if (angle_type[atom2][i] > 0)
		energy[m] = angle->single(angle_type[atom2][i],atom1,atom2,atom3);
		else energy[m] = 0.0;
		ebuf[n] = angle->single(angle_type[atom2][i],atom1,atom2,atom3);
		else ebuf[n] = 0.0;
		}
		n += nvalues;
		}

		m++;
		@@ -200,46 +190,17 @@ int ComputeAngleLocal::compute_angles(int flag)

		/* ---------------------------------------------------------------------- */

		void ComputeAngleLocal::pack_theta(int n)
		{
		for (int m = 0; m < ncount; m++) {
		buf[n] = theta[m];
		n += nvalues;
		}
		}

		/* ---------------------------------------------------------------------- */

		void ComputeAngleLocal::pack_energy(int n)
		{
		for (int m = 0; m < ncount; m++) {
		buf[n] = energy[m];
		n += nvalues;
		}
		}

		/* ---------------------------------------------------------------------- */

		void ComputeAngleLocal::reallocate(int n)
		{
		// grow vector or array and indices array

		while (nmax < n) nmax += DELTA;

		if (tflag) {
		memory->sfree(theta);
		theta = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:theta");
		}
		if (eflag) {
		memory->sfree(energy);
		energy = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:energy");
		}

		if (nvalues == 1) {
		if (tflag) vector_local = theta;
		if (eflag) vector_local = energy;
		memory->sfree(vector);
		vector = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:vector");
		vector_local = vector;
		} else {
		memory->destroy_2d_double_array(array);
		array = memory->create_2d_double_array(nmax,nvalues,
		@@ -254,9 +215,6 @@ void ComputeAngleLocal::reallocate(int n)

		double ComputeAngleLocal::memory_usage()
		{
		double bytes = 0.0;
		if (tflag) bytes += nmax * sizeof(double);
		if (eflag) bytes += nmax * sizeof(double);
		if (nvalues > 1) bytes += nmaxnvalues sizeof(double);
		double bytes = nmaxnvalues sizeof(double);
		return bytes;
		}

src/compute_angle_local.h

+1 −10

Original line number	Diff line number	Diff line
		@@ -28,23 +28,14 @@ class ComputeAngleLocal : public Compute {

		private:
		int nvalues,tflag,eflag;
		int *which;
		int ncount;

		int nmax;
		double *theta;
		double *energy;
		double *vector;
		double **array;
		double *buf;

		int compute_angles(int);
		void reallocate(int);

		typedef void (ComputeAngleLocal::*FnPtrPack)(int);
		FnPtrPack *pack_choice; // ptrs to pack functions

		void pack_theta(int);
		void pack_energy(int);
		};

		}

src/compute_bond_local.cpp

+30 −72

Original line number	Diff line number	Diff line
		@@ -25,8 +25,6 @@

		using namespace LAMMPS_NS;

		enum{DISTANCE,ENERGY};

		#define DELTA 10000

		/* ---------------------------------------------------------------------- */
		@@ -44,28 +42,19 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS lmp, int narg, char *arg) :
		if (nvalues == 1) size_local_cols = 0;
		else size_local_cols = nvalues;

		which = new int[nvalues];
		pack_choice = new FnPtrPack[nvalues];

		dflag = eflag = 0;
		dflag = eflag = -1;
		nvalues = 0;

		int i;
		for (int iarg = 3; iarg < narg; iarg++) {
		i = iarg-3;

		if (strcmp(arg[iarg],"distance") == 0) {
		dflag = 1;
		which[i] = DISTANCE;
		pack_choice[i] = &ComputeBondLocal::pack_distance;
		} else if (strcmp(arg[iarg],"energy") == 0) {
		eflag = 1;
		which[i] = ENERGY;
		pack_choice[i] = &ComputeBondLocal::pack_energy;
		} else error->all("Invalid keyword in compute bond/local command");
		if (strcmp(arg[iarg],"dist") == 0) dflag = nvalues++;
		else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
		else error->all("Invalid keyword in compute bond/local command");
		}

		nmax = 0;
		distance = energy = NULL;
		vector = NULL;
		array = NULL;
		}

		@@ -73,10 +62,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS lmp, int narg, char *arg) :

		ComputeBondLocal::~ComputeBondLocal()
		{
		delete [] which;
		delete [] pack_choice;
		memory->sfree(distance);
		memory->sfree(energy);
		memory->sfree(vector);
		memory->destroy_2d_double_array(array);
		}

		@@ -106,14 +92,6 @@ void ComputeBondLocal::compute_local()
		if (ncount > nmax) reallocate(ncount);
		size_local_rows = ncount;
		ncount = compute_bonds(1);

		// fill array with distance/energy values

		if (nvalues > 1) {
		if (array) buf = array[0];
		for (int n = 0; n < nvalues; n++)
		(this->*pack_choice[n])(n);
		}
		}

		/* ----------------------------------------------------------------------
		@@ -129,8 +107,9 @@ void ComputeBondLocal::compute_local()

		int ComputeBondLocal::compute_bonds(int flag)
		{
		int i,atom1,atom2;
		int i,m,n,atom1,atom2;
		double delx,dely,delz,rsq;
		double dbuf,ebuf;

		double **x = atom->x;
		int *num_bond = atom->num_bond;
		@@ -141,9 +120,19 @@ int ComputeBondLocal::compute_bonds(int flag)
		int nlocal = atom->nlocal;
		int newton_bond = force->newton_bond;

		if (flag) {
		if (nvalues == 1) {
		if (dflag >= 0) dbuf = vector;
		if (eflag >= 0) ebuf = vector;
		} else {
		if (dflag >= 0) dbuf = &array[0][dflag];
		if (eflag >= 0) ebuf = &array[0][eflag];
		}
		}

		Bond *bond = force->bond;

		int m = 0;
		m = n = 0;
		for (atom1 = 0; atom1 < nlocal; atom1++) {
		if (!(mask[atom1] & groupbit)) continue;
		for (i = 0; i < num_bond[atom1]; i++) {
		@@ -158,12 +147,13 @@ int ComputeBondLocal::compute_bonds(int flag)
		delz = x[atom1][2] - x[atom2][2];
		domain->minimum_image(delx,dely,delz);
		rsq = delxdelx + delydely + delz*delz;
		if (dflag) distance[m] = sqrt(rsq);
		if (eflag) {
		if (dflag >= 0) dbuf[n] = sqrt(rsq);
		if (eflag >= 0) {
		if (bond_type[atom1][i] > 0)
		energy[m] = bond->single(bond_type[atom1][i],rsq,atom1,atom2);
		else energy[m] = 0.0;
		ebuf[n] = bond->single(bond_type[atom1][i],rsq,atom1,atom2);
		else ebuf[n] = 0.0;
		}
		n += nvalues;
		}

		m++;
		@@ -175,46 +165,17 @@ int ComputeBondLocal::compute_bonds(int flag)

		/* ---------------------------------------------------------------------- */

		void ComputeBondLocal::pack_distance(int n)
		{
		for (int m = 0; m < ncount; m++) {
		buf[n] = distance[m];
		n += nvalues;
		}
		}

		/* ---------------------------------------------------------------------- */

		void ComputeBondLocal::pack_energy(int n)
		{
		for (int m = 0; m < ncount; m++) {
		buf[n] = energy[m];
		n += nvalues;
		}
		}

		/* ---------------------------------------------------------------------- */

		void ComputeBondLocal::reallocate(int n)
		{
		// grow vector or array and indices array

		while (nmax < n) nmax += DELTA;

		if (dflag) {
		memory->sfree(distance);
		distance = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:distance");
		}
		if (eflag) {
		memory->sfree(energy);
		energy = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:energy");
		}

		if (nvalues == 1) {
		if (dflag) vector_local = distance;
		if (eflag) vector_local = energy;
		memory->sfree(vector);
		vector = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:vector");
		vector_local = vector;
		} else {
		memory->destroy_2d_double_array(array);
		array = memory->create_2d_double_array(nmax,nvalues,
		@@ -229,9 +190,6 @@ void ComputeBondLocal::reallocate(int n)

		double ComputeBondLocal::memory_usage()
		{
		double bytes = 0.0;
		if (dflag) bytes += nmax * sizeof(double);
		if (eflag) bytes += nmax * sizeof(double);
		if (nvalues > 1) bytes += nmaxnvalues sizeof(double);
		double bytes = nmaxnvalues sizeof(double);
		return bytes;
		}

src/compute_bond_local.h

+1 −9

Original line number	Diff line number	Diff line
		@@ -28,23 +28,15 @@ class ComputeBondLocal : public Compute {

		private:
		int nvalues,dflag,eflag;
		int *which;
		int ncount;

		int nmax;
		double *distance;
		double *energy;
		double *vector;
		double **array;
		double *buf;

		int compute_bonds(int);
		void reallocate(int);

		typedef void (ComputeBondLocal::*FnPtrPack)(int);
		FnPtrPack *pack_choice; // ptrs to pack functions

		void pack_distance(int);
		void pack_energy(int);
		};

		}

src/compute_dihedral_local.cpp

0 → 100644

+256 −0

Original line number	Diff line number	Diff line
		/* ----------------------------------------------------------------------
		LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
		http://lammps.sandia.gov, Sandia National Laboratories
		Steve Plimpton, sjplimp@sandia.gov

		Copyright (2003) Sandia Corporation. Under the terms of Contract
		DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
		certain rights in this software. This software is distributed under
		the GNU General Public License.

		See the README file in the top-level LAMMPS directory.
		------------------------------------------------------------------------- */

		#include "math.h"
		#include "string.h"
		#include "compute_dihedral_local.h"
		#include "atom.h"
		#include "atom_vec.h"
		#include "update.h"
		#include "domain.h"
		#include "force.h"
		#include "dihedral.h"
		#include "memory.h"
		#include "error.h"

		using namespace LAMMPS_NS;

		#define DELTA 10000

		#define MIN(A,B) ((A) < (B)) ? (A) : (B)
		#define MAX(A,B) ((A) > (B)) ? (A) : (B)
		#define SMALL 0.001

		/* ---------------------------------------------------------------------- */

		ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS lmp, int narg, char *arg) :
		Compute(lmp, narg, arg)
		{
		if (narg < 4) error->all("Illegal compute dihedral/local command");

		if (atom->avec->dihedrals_allow == 0)
		error->all("Compute dihedral/local used when dihedrals are not allowed");

		local_flag = 1;
		nvalues = narg - 3;
		if (nvalues == 1) size_local_cols = 0;
		else size_local_cols = nvalues;

		pflag = -1;
		nvalues = 0;

		int i;
		for (int iarg = 3; iarg < narg; iarg++) {
		i = iarg-3;
		if (strcmp(arg[iarg],"phi") == 0) pflag = nvalues++;
		else error->all("Invalid keyword in compute dihedral/local command");
		}

		nmax = 0;
		vector = NULL;
		array = NULL;
		}

		/* ---------------------------------------------------------------------- */

		ComputeDihedralLocal::~ComputeDihedralLocal()
		{
		memory->sfree(vector);
		memory->destroy_2d_double_array(array);
		}

		/* ---------------------------------------------------------------------- */

		void ComputeDihedralLocal::init()
		{
		if (force->dihedral == NULL)
		error->all("No dihedral style is defined for compute dihedral/local");

		// do initial memory allocation so that memory_usage() is correct

		ncount = compute_dihedrals(0);
		if (ncount > nmax) reallocate(ncount);
		size_local_rows = ncount;
		}

		/* ---------------------------------------------------------------------- */

		void ComputeDihedralLocal::compute_local()
		{
		invoked_local = update->ntimestep;

		// count local entries and compute dihedral info

		ncount = compute_dihedrals(0);
		if (ncount > nmax) reallocate(ncount);
		size_local_rows = ncount;
		ncount = compute_dihedrals(1);
		}

		/* ----------------------------------------------------------------------
		count dihedrals on this proc
		only count if 2nd atom is the one storing the dihedral
		all atoms in interaction must be in group
		all atoms in interaction must be known to proc
		if flag is set, compute requested info about dihedral
		------------------------------------------------------------------------- */

		int ComputeDihedralLocal::compute_dihedrals(int flag)
		{
		int i,m,n,atom1,atom2,atom3,atom4;
		double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
		double sb1,sb2,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2;
		double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2;
		double c2mag,sin2,sc1,sc2,s1,s2,s12,c;
		double *pbuf;

		double **x = atom->x;
		int *num_dihedral = atom->num_dihedral;
		int **dihedral_atom1 = atom->dihedral_atom1;
		int **dihedral_atom2 = atom->dihedral_atom2;
		int **dihedral_atom3 = atom->dihedral_atom3;
		int **dihedral_atom4 = atom->dihedral_atom4;
		int *tag = atom->tag;
		int *mask = atom->mask;
		int nlocal = atom->nlocal;

		if (flag) {
		if (nvalues == 1) {
		if (pflag >= 0) pbuf = vector;
		} else {
		if (pflag >= 0) pbuf = &array[0][pflag];
		}
		}

		double PI = 4.0*atan(1.0);

		m = n = 0;
		for (atom2 = 0; atom2 < nlocal; atom2++) {
		if (!(mask[atom2] & groupbit)) continue;
		for (i = 0; i < num_dihedral[atom2]; i++) {
		if (tag[atom2] != dihedral_atom2[atom2][i]) continue;
		atom1 = atom->map(dihedral_atom1[atom2][i]);
		if (atom1 < 0 \|\| !(mask[atom1] & groupbit)) continue;
		atom3 = atom->map(dihedral_atom3[atom2][i]);
		if (atom3 < 0 \|\| !(mask[atom3] & groupbit)) continue;
		atom4 = atom->map(dihedral_atom4[atom2][i]);
		if (atom4 < 0 \|\| !(mask[atom4] & groupbit)) continue;

		if (flag) {

		// phi calculation from dihedral style OPLS

		if (pflag >= 0) {
		vb1x = x[atom1][0] - x[atom2][0];
		vb1y = x[atom1][1] - x[atom2][1];
		vb1z = x[atom1][2] - x[atom2][2];
		domain->minimum_image(vb1x,vb1y,vb1z);

		vb2x = x[atom3][0] - x[atom2][0];
		vb2y = x[atom3][1] - x[atom2][1];
		vb2z = x[atom3][2] - x[atom2][2];
		domain->minimum_image(vb2x,vb2y,vb2z);

		vb2xm = -vb2x;
		vb2ym = -vb2y;
		vb2zm = -vb2z;
		domain->minimum_image(vb2xm,vb2ym,vb2zm);

		vb3x = x[atom4][0] - x[atom3][0];
		vb3y = x[atom4][1] - x[atom3][1];
		vb3z = x[atom4][2] - x[atom3][2];
		domain->minimum_image(vb3x,vb3y,vb3z);

		sb1 = 1.0 / (vb1xvb1x + vb1yvb1y + vb1z*vb1z);
		sb2 = 1.0 / (vb2xvb2x + vb2yvb2y + vb2z*vb2z);
		sb3 = 1.0 / (vb3xvb3x + vb3yvb3y + vb3z*vb3z);

		rb1 = sqrt(sb1);
		rb3 = sqrt(sb3);

		c0 = (vb1xvb3x + vb1yvb3y + vb1zvb3z) rb1*rb3;

		b1mag2 = vb1xvb1x + vb1yvb1y + vb1z*vb1z;
		b1mag = sqrt(b1mag2);
		b2mag2 = vb2xvb2x + vb2yvb2y + vb2z*vb2z;
		b2mag = sqrt(b2mag2);
		b3mag2 = vb3xvb3x + vb3yvb3y + vb3z*vb3z;
		b3mag = sqrt(b3mag2);

		ctmp = vb1xvb2x + vb1yvb2y + vb1z*vb2z;
		r12c1 = 1.0 / (b1mag*b2mag);
		c1mag = ctmp * r12c1;

		ctmp = vb2xmvb3x + vb2ymvb3y + vb2zm*vb3z;
		r12c2 = 1.0 / (b2mag*b3mag);
		c2mag = ctmp * r12c2;

		sin2 = MAX(1.0 - c1mag*c1mag,0.0);
		sc1 = sqrt(sin2);
		if (sc1 < SMALL) sc1 = SMALL;
		sc1 = 1.0/sc1;

		sin2 = MAX(1.0 - c2mag*c2mag,0.0);
		sc2 = sqrt(sin2);
		if (sc2 < SMALL) sc2 = SMALL;
		sc2 = 1.0/sc2;

		s1 = sc1 * sc1;
		s2 = sc2 * sc2;
		s12 = sc1 * sc2;
		c = (c0 + c1magc2mag) s12;

		if (c > 1.0) c = 1.0;
		if (c < -1.0) c = -1.0;
		pbuf[n] = 180.0*acos(c)/PI;
		}
		n += nvalues;
		}

		m++;
		}
		}

		return m;
		}

		/* ---------------------------------------------------------------------- */

		void ComputeDihedralLocal::reallocate(int n)
		{
		// grow vector or array and indices array

		while (nmax < n) nmax += DELTA;

		if (nvalues == 1) {
		memory->sfree(vector);
		vector = (double ) memory->smalloc(nmaxsizeof(double),
		"bond/local:vector");
		vector_local = vector;
		} else {
		memory->destroy_2d_double_array(array);
		array = memory->create_2d_double_array(nmax,nvalues,
		"bond/local:array");
		array_local = array;
		}
		}

		/* ----------------------------------------------------------------------
		memory usage of local data
		------------------------------------------------------------------------- */

		double ComputeDihedralLocal::memory_usage()
		{
		double bytes = nmaxnvalues sizeof(double);
		return bytes;
		}

Admin message