</div>
<div class="method-description">
<p>
-Find atoms that are within the threshold distance from the given atom. If limit is a positive number, the threshold distance is the sum of the vdw radii times limit. If limit is a negative number, its absolute value is used for the threshold distance in angstrom. If limit is not given, a default value of 1.2 is used. The number of the newly created bonds is returned. An array of atom indices is returned. If no atoms are found, an empty array is returned.
+Find atoms that are within the threshold distance from the given atom. If limit is a positive number, the threshold distance is the sum of the vdw radii times limit. If limit is a negative number, its absolute value is used for the threshold distance in angstrom. If limit is not given, a default value of 1.2 is used.
+</p>
+<p>
+An array of atom indices is returned. If no atoms are found, an empty array is returned.
+</p>
+<p>
+The atom argument can also be a <a href="Vector3D.html">Vector3D</a> (or any other object that can be converted to a <a href="Vector3D.html">Vector3D</a>), representing cartesian coordinates.
</p>
</div>
</div>
#include "MDGraphite.h"
#include "MDPressure.h"
#include "MDEwald.h"
+#include "MDSurface.h"
+#include "MDForce.h"
#include <stdlib.h>
#include <string.h>
void md_init_for_positions(MDArena *arena);
const char *md_prepare(MDArena *arena, int check_only);
+int md_set_external_forces(MDArena *arena, int nexforces, const Vector *exforces);
int md_check_abnormal_bond(MDArena *arena, Molecule *mol, int idx);
int md_check_abnormal_angle(MDArena *arena, Molecule *mol, int idx);
int md_check_abnormal_dihedral(MDArena *arena, Molecule *mol, int idx);
int md_output_results(MDArena *arena);
int md_copy_coordinates_from_internal(MDArena *arena);
+int md_copy_coordinates_to_internal(MDArena *arena);
int md_is_running(MDArena *arena);
void calc_force(MDArena *arena);
void calc_pair_interaction(MDArena *arena, const MDGroupFlags *group_flags_1, const MDGroupFlags *group_flags_2);
-
+void calc_surface_force(MDArena *arena);
+
#ifdef __cplusplus
}
#endif
*/
#include "MDPressure.h"
+#include "MDForce.h"
#include <stdlib.h>
#include <math.h>
the rest of MolLib.h. */
#include "MainView.h"
#include "MolLib.h"
+#include "Molby_extern.h"
#include "MD/MDCore.h"
#include "MD/MDGraphite.h"
case 'i': lastIval = argp->u.ival; fprintf(fp, " %d", lastIval); break;
case 'd': fprintf(fp, " %g", argp->u.dval); break;
case 's':
- case 'C': fprintf(fp, " %.*s", argp->u.arval.nitems, argp->u.arval.ptr); break;
+ case 'C': fprintf(fp, " %.*s", argp->u.arval.nitems, (char *)argp->u.arval.ptr); break;
case 'v': case 't': case 'u': case 'I': case 'D': case 'V': case 'T': case 'U': {
char *pp = (char *)argp->u.arval.ptr;
int n = argp->u.arval.nitems;
MoleculeCallback_displayName(argp->u.mval, buf, sizeof buf);
if (buf[0] == 0) {
/* No associated document */
- snprintf(buf, sizeof buf, "#<Molecule:0x%lx>", argp->u.mval);
+ snprintf(buf, sizeof buf, "#<Molecule:0x%lx>", (long unsigned)(void *)argp->u.mval);
}
fprintf(fp, "%s", buf);
break;
j = ParameterReadFromString(par, buf, &bufp, fname, lineNumber, 0);
if (j < 0) {
s_append_asprintf(errbuf, "%s", bufp);
+ free(bufp);
goto err_exit;
}
i += j;
}
if (bufp != NULL) {
- MyAppCallback_setConsoleColor(1);
- MyAppCallback_showScriptMessage("%s", bufp);
- MyAppCallback_setConsoleColor(0);
+ s_append_asprintf(errbuf, "%s", bufp);
free(bufp);
}
continue;
/* Subroutine for MoleculeGuessBonds. It can be also used independently, but make sure that *outNbonds/*outBonds
_correctly_ represents an array of two integers (as in mp->nbonds/mp->bonds). */
-/* Find atoms within the given "distance" from the given atom. */
+/* Find atoms within the given "distance" from the given position. */
/* If limit is negative, its absolute value denotes the threshold distance in angstrom; otherwise,
- the threshold distance is given by the sum of van der Waals radii times limit. */
-/* If triangle is non-zero, then only atoms with lower indexes than index are looked for. */
+ the threshold distance is given by the sum of van der Waals radii times limit, and radius is
+ the van der Waals radius of the atom at the given position. */
+/* Index is the atom index of the given atom; it is only used in returning the "bond" array
+ to the caller. If index is negative, then (-index) is the real atom index, and
+ only atoms with lower indices than (-index) are looked for. */
int
-MoleculeFindCloseAtoms(Molecule *mp, Int index, Double limit, Int *outNbonds, Int **outBonds, Int triangle)
-{
- Int n1, n2, j, nlim, newbond[2];
- Double a1, a2, alim;
- Vector dr, r1, r2;
- Atom *ap = ATOM_AT_INDEX(mp->atoms, index);
- n1 = ap->atomicNumber;
- if (n1 >= 0 && n1 < gCountElementParameters)
- a1 = gElementParameters[n1].radius;
- else a1 = gElementParameters[6].radius;
- r1 = ap->r;
- nlim = (triangle ? index : mp->natoms);
+MoleculeFindCloseAtoms(Molecule *mp, const Vector *vp, Double radius, Double limit, Int *outNbonds, Int **outBonds, Int index)
+{
+ Int n2, j, nlim, newbond[2];
+ Double a2, alim;
+ Vector dr, r2;
+ if (index < 0) {
+ nlim = index = -index;
+ } else {
+ nlim = mp->natoms;
+ }
for (j = 0; j < nlim; j++) {
Atom *bp = ATOM_AT_INDEX(mp->atoms, j);
if (index == j)
a2 = gElementParameters[n2].radius;
else a2 = gElementParameters[6].radius;
r2 = bp->r;
- VecSub(dr, r1, r2);
+ VecSub(dr, *vp, r2);
if (limit < 0)
alim = -limit;
else
- alim = limit * (a1 + a2);
+ alim = limit * (radius + a2);
if (VecLength2(dr) < alim * alim) {
newbond[0] = index;
newbond[1] = j;
MoleculeGuessBonds(Molecule *mp, Double limit, Int *outNbonds, Int **outBonds)
{
Int nbonds, *bonds, i, newbond[2];
-/* int i, j, n1, n2;
- Atom *ap, *bp;
- Vector r1, r2, dr;
- Double a1, a2, alim;
- Int newbond[2];
- ElementPar *p = gElementParameters; */
+ Atom *ap;
nbonds = 0;
bonds = NULL;
- for (i = 0; i < mp->natoms; i++) {
- MoleculeFindCloseAtoms(mp, i, limit, &nbonds, &bonds, 1);
- /*
- ap = ATOM_AT_INDEX(mp->atoms, i);
- n1 = ap->atomicNumber;
- if (n1 >= 0 && n1 < gCountElementParameters)
- a1 = p[n1].radius;
- else a1 = p[6].radius;
- r1 = ap->r;
- for (j = 0; j < i; j++) {
- bp = ATOM_AT_INDEX(mp->atoms, j);
- n2 = bp->atomicNumber;
- if (n2 >= 0 && n2 < gCountElementParameters)
- a2 = p[n2].radius;
- else a2 = p[6].radius;
- r2 = bp->r;
- VecSub(dr, r1, r2);
- if (limit < 0)
- alim = -limit;
- else
- alim = limit * (a1 + a2);
- if (VecLength2(dr) < alim * alim) {
- newbond[0] = i;
- newbond[1] = j;
- AssignArray(&bonds, &nbonds, sizeof(Int) * 2, nbonds, newbond);
- }
- }
- */
+ for (i = 1, ap = ATOM_NEXT(mp->atoms); i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ Vector r = ap->r;
+ Int an = ap->atomicNumber;
+ Double rad;
+ if (an >= 0 && an < gCountElementParameters)
+ rad = gElementParameters[an].radius;
+ else rad = gElementParameters[6].radius;
+ MoleculeFindCloseAtoms(mp, &r, rad, limit, &nbonds, &bonds, -i);
}
if (nbonds > 0) {
newbond[0] = kInvalidIndex;
void MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double x33, Double x12, Double x13, Double x23, const Double *sigmap);
void MoleculeSetAnisoBySymop(Molecule *mp, int idx);
int MoleculeSetPeriodicBox(Molecule *mp, const Vector *ax, const Vector *ay, const Vector *az, const Vector *ao, const char *periodic, int convertCoordinates);
+int MoleculeCalculateCellFromAxes(XtalCell *cp, int calc_abc);
int MoleculeReadCoordinatesFromFile(Molecule *mp, const char *fname, const char *ftype, char **errbuf);
int MoleculeReadCoordinatesFromPdbFile(Molecule *mp, const char *fname, char **errbuf);
int MoleculeReadCoordinatesFromDcdFile(Molecule *mp, const char *fname, char **errbuf);
+int MoleculeLoadGamessDatFile(Molecule *mol, const char *fname, char **errbuf);
+
int MoleculeReadExtendedInfo(Molecule *mp, const char *fname, char **errbuf);
int MoleculeWriteExtendedInfo(Molecule *mp, const char *fname, char **errbuf);
int MoleculeWriteToPdbFile(Molecule *mp, const char *fname, char **errbuf);
int MoleculeWriteToDcdFile(Molecule *mp, const char *fname, char **errbuf);
int MoleculeWriteToTepFile(Molecule *mp, const char *fname, char **errbuf);
+int MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf);
void MoleculeDump(Molecule *mol);
int MoleculePrepareMDArena(Molecule *mol, int check_only, char **retmsg);
int MoleculeAnalyzeAtomName(const char *s, char *resName, int *resSeq, char *atomName);
int MoleculeAtomIndexFromString(Molecule *mp, const char *s);
-int MoleculeFindCloseAtoms(Molecule *mp, Int index, Double limit, Int *outNbonds, Int **outBonds, Int triangle);
+int MoleculeFindCloseAtoms(Molecule *mp, const Vector *vp, Double radius, Double limit, Int *outNbonds, Int **outBonds, Int triangle);
int MoleculeGuessBonds(Molecule *mp, Double limit, Int *outNbonds, Int **outBonds);
int MoleculeRebuildTablesFromConnects(Molecule *mp);
int MoleculeAreAtomsConnected(Molecule *mol, int idx1, int idx2);
int ParameterGlobalParIndexForSrcIndex(int src);
int ParameterCommentIndexForGlobalFileName(const char *p);
+int ParameterCompare(const UnionPar *up1, const UnionPar *up2, int type);
int ParameterCommentIndex(const char *comment);
const char *ParameterGetComment(int idx);
extern VALUE ValueFromIntGroup(IntGroup *ig);
extern VALUE ValueFromVector(const Vector *vp);
extern VALUE ValueFromTransform(Transform *tp);
+extern VALUE ValueFromMoleculeWithParameterTypeAndIndex(Molecule *mol, int type, int idx1);
extern void VectorFromValue(VALUE val, Vector *vp);
extern void TransformFromValue(VALUE val, Transform *tp);
/*
* call-seq:
- * find_close_atoms(atom, limit = 1.2) -> array of Integers (atom indices)
+ * find_close_atoms(atom, limit = 1.2, radius = 0.77) -> array of Integers (atom indices)
*
* Find atoms that are within the threshold distance from the given atom.
+ * (The atom argument can also be a vector, representing a cartesian coordinate. In that case, the van der Waals of the atom can also be specified.)
* If limit is a positive number, the threshold distance is the sum of the vdw radii times limit.
* If limit is a negative number, its absolute value is used for the threshold distance in angstrom.
* If limit is not given, a default value of 1.2 is used.
s_Molecule_FindCloseAtoms(int argc, VALUE *argv, VALUE self)
{
Molecule *mol;
- VALUE aval, limval;
- double limit;
- Int n1, nbonds, *bonds;
+ VALUE aval, limval, radval;
+ double limit, radius;
+ Int n1, nbonds, *bonds, an;
+ Vector v;
Data_Get_Struct(self, Molecule, mol);
- rb_scan_args(argc, argv, "11", &aval, &limval);
- n1 = s_Molecule_AtomIndexFromValue(mol, aval);
+ rb_scan_args(argc, argv, "12", &aval, &limval, &radval);
+ if (rb_obj_is_kind_of(aval, rb_cVector3D) || rb_obj_is_kind_of(aval, rb_cLAMatrix) || rb_obj_is_kind_of(aval, rb_mEnumerable)) {
+ VectorFromValue(aval, &v);
+ if (radval == Qnil)
+ radius = gElementParameters[6].radius;
+ else
+ radius = NUM2DBL(rb_Float(radval));
+ n1 = 0;
+ } else {
+ n1 = s_Molecule_AtomIndexFromValue(mol, aval);
+ v = ATOM_AT_INDEX(mol->atoms, n1)->r;
+ an = ATOM_AT_INDEX(mol->atoms, n1)->atomicNumber;
+ if (an >= 0 && an < gCountElementParameters)
+ radius = gElementParameters[an].radius;
+ else radius = gElementParameters[6].radius;
+ }
if (limval == Qnil)
limit = 1.2;
else
limit = NUM2DBL(rb_Float(limval));
nbonds = 0; /* This initialization is necessary: see comments in MoleculeFindCloseAtoms() */
bonds = NULL;
- MoleculeFindCloseAtoms(mol, n1, limit, &nbonds, &bonds, 0);
+ MoleculeFindCloseAtoms(mol, &v, radius, limit, &nbonds, &bonds, n1);
aval = rb_ary_new();
if (nbonds > 0) {
for (n1 = 0; n1 < nbonds; n1++)
static VALUE
s_MDArena_testPME(int argc, VALUE *argv, VALUE self)
{
- extern pme_test(MDArena *);
+ extern int pme_test(MDArena *);
MDArena *arena;
if (rb_obj_is_kind_of(self, rb_cMDArena)) {
Data_Get_Struct(self, MDArena, arena);
/* Avoid this "dummy" top-level window to appear in the window menu.
It should not happen because MyApp::OnActivate() tries to hide this window,
but this is still here just in case. */
- OSStatus sts;
+// OSStatus sts;
// sts = ChangeWindowAttributes((WindowRef)m_macWindow, 0, kWindowInWindowMenuAttribute);
/* printf("m_macWindow = %p, status = %d\n", m_macWindow, (int)sts); */
#endif
MyGLCanvas::OnSize(wxSizeEvent &event)
{
// this is also necessary to update the context on some platforms
- wxGLCanvas::OnSize(event);
+// wxGLCanvas::OnSize(event);
// set GL viewport (not called by wxGLCanvas::OnSize on all platforms...)
int w, h;
if (m_exit_func)
(*m_exit_func)(m_argptr, m_argnum);
Exit(code);
+ return 0;
}
OSDN Git Service
