/* Check the bonded atoms and append to results if not already present */
/* results[] is terminated by -1, hence must be at least (natom+1) size */
static int
-s_check_bonded(Atom *ap, Int *results)
+s_check_bonded(Molecule *mol, Int idx, Int *results)
{
Int i, n, *ip;
const Int *cp;
+ Atom *ap = ATOM_AT_INDEX(mol->atoms, idx);
cp = AtomConnectData(&ap->connect);
for (i = 0; i < ap->connect.count; i++, cp++) {
n = *cp;
*ip = -1;
}
}
+ if (mol->piconnects != NULL) {
+ for (i = mol->piconnects[idx]; i < mol->piconnects[idx + 1]; i++) {
+ n = mol->piconnects[i];
+ for (ip = results; *ip >= 0; ip++) {
+ if (n == *ip)
+ break;
+ }
+ if (*ip < 0) {
+ *ip++ = n;
+ *ip = -1;
+ }
+ }
+ }
for (ip = results; *ip >= 0; ip++)
;
return ip - results;
{
Int *results;
Int natoms = arena->mol->natoms;
- Atom *atoms = arena->mol->atoms;
MDExclusion *exinfo;
int next_index, i, j;
results[1] = -1;
exinfo[i].index1 = 1;
/* 1-2 exclusion (directly bonded) */
- exinfo[i].index2 = s_check_bonded(&atoms[i], results);
+ exinfo[i].index2 = s_check_bonded(arena->mol, i, results);
n = exinfo[i].index2;
/* 1-3 exclusion: atoms bonded to 1-2 exclusions */
for (j = exinfo[i].index1; j < exinfo[i].index2; j++)
- n = s_check_bonded(&atoms[results[j]], results);
+ n = s_check_bonded(arena->mol, results[j], results);
exinfo[i].index3 = n;
/* 1-4 exclusion: atoms bonded to 1-3 exclusions */
for (j = exinfo[i].index2; j < exinfo[i].index3; j++)
- n = s_check_bonded(&atoms[results[j]], results);
+ n = s_check_bonded(arena->mol, results[j], results);
AssignArray(&arena->exlist, &arena->nexlist, sizeof(Int), next_index + n, NULL);
memcpy(arena->exlist + next_index, results, n * sizeof(Int));
exinfo[i].index0 += next_index;
ptype = -1;
break;
} else
- types[j] = mol->piatoms[j].type;
+ types[j] = mol->piatoms[ii].type;
} else if (ii < mol->natoms) {
idx[j] = ii;
types[j] = ATOM_AT_INDEX(mol->atoms, ii)->type;
NewArray(&(mp2->pibonds), &(mp2->npibonds), sizeof(Int) * 4, mp->npibonds);
memmove(mp2->pibonds, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
}
+ if (mp->npiconnects > 0) {
+ NewArray(&(mp2->piconnects), &(mp2->npiconnects), sizeof(Int), mp->npiconnects);
+ memmove(mp2->piconnects, mp->piconnects, sizeof(Int) * mp->npiconnects);
+ }
/* mp2->useFlexibleCell = mp->useFlexibleCell; */
if (mp->nframe_cells > 0) {
mp->pibonds = NULL;
mp->npibonds = 0;
}
+ if (mp->piconnects != NULL) {
+ free(mp->piconnects);
+ mp->piconnects = NULL;
+ mp->npiconnects = 0;
+ }
if (mp->selection != NULL) {
IntGroupRelease(mp->selection);
mp->selection = NULL;
return -2;
/* Create atoms, with avoiding duplicates */
+ MoleculeInvalidatePiConnectionTable(mp);
n0 = n1 = mp->natoms;
table = (int *)malloc(sizeof(int) * n0);
if (table == NULL)
if (mp->impropers[i] >= pos)
mp->impropers[i]++;
}
+ MoleculeInvalidatePiConnectionTable(mp);
}
mp->nframes = -1; /* Should be recalculated later */
MoleculeIncrementModifyCount(mp);
return 1; /* Bad parameter */
__MoleculeLock(dst);
+ MoleculeInvalidatePiConnectionTable(dst);
nsrc = src->natoms;
ndst = dst->natoms;
if (resSeqOffset < 0)
return 1; /* Bad parameter */
__MoleculeLock(src);
+ MoleculeInvalidatePiConnectionTable(src);
nsrc = src->natoms;
nsrcnew = nsrc - ndst;
free(newAtoms);
free(old2new);
free(apArray);
+
+ MoleculeInvalidatePiConnectionTable(mp);
}
for (i = 0; i < mp->natoms; i++)
memmove(ATOM_AT_INDEX(mp->atoms, old2new[i]), ATOM_AT_INDEX(saveAtoms, i), gSizeOfAtomRecord);
retval = 0;
+
+ MoleculeInvalidatePiConnectionTable(mp);
MoleculeIncrementModifyCount(mp);
mp->needsMDRebuild = 1;
sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
{
Atom *ap;
- Int i, *cp;
+ Int i, *cp, idx2;
if (exatoms != NULL && IntGroupLookup(exatoms, idx, NULL))
return;
IntGroupAdd(result, idx, 1);
ap = ATOM_AT_INDEX(mp->atoms, idx);
cp = AtomConnectData(&ap->connect);
for (i = 0; i < ap->connect.count; i++) {
- int idx2 = cp[i];
+ idx2 = cp[i];
if (IntGroupLookup(result, idx2, NULL))
continue;
sMoleculeFragmentSub(mp, idx2, result, exatoms);
}
+ if (mp->piconnects != NULL) {
+ for (i = mp->piconnects[idx]; i < mp->piconnects[idx + 1]; i++) {
+ idx2 = mp->piconnects[i];
+ if (IntGroupLookup(result, idx2, NULL))
+ continue;
+ sMoleculeFragmentSub(mp, idx2, result, exatoms);
+ }
+ }
}
/* The molecular fragment (= interconnected atoms) containing the atom n1 and
if (mp == NULL || mp->natoms == 0 || n1 < 0 || n1 >= mp->natoms)
return NULL;
result = IntGroupNew();
+ MoleculeValidatePiConnectionTable(mp);
sMoleculeFragmentSub(mp, n1, result, exatoms);
return result;
}
for (i = 0; i < argc; i++)
IntGroupAdd(exatoms, argv[i], 1);
result = IntGroupNew();
+ MoleculeValidatePiConnectionTable(mp);
sMoleculeFragmentSub(mp, n1, result, exatoms);
IntGroupRelease(exatoms);
return result;
return NULL;
IntGroupIteratorInit(inatoms, &iter);
result = IntGroupNew();
+ MoleculeValidatePiConnectionTable(mp);
while ((i = IntGroupIteratorNext(&iter)) >= 0) {
sMoleculeFragmentSub(mp, i, result, exatoms);
}
int
MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
{
- Int i, i1, i2, j, k, bond_count, nval1, nval2, *cp;
+ Int i, i1, i2, j, k, k2, bond_count, nval1, nval2, *cp;
Atom *ap;
if (mp == NULL || mp->natoms == 0 || group == NULL)
return 0; /* Invalid arguments */
+ MoleculeValidatePiConnectionTable(mp);
bond_count = 0;
for (i = 0; (i1 = IntGroupGetStartPoint(group, i)) >= 0; i++) {
i2 = IntGroupGetEndPoint(group, i);
return 0; /* Too many bonds */
}
}
+ if (mp->piconnects != NULL) {
+ for (k = mp->piconnects[j]; k < mp->piconnects[j + 1]; k++) {
+ k2 = mp->piconnects[k];
+ if (IntGroupLookup(group, k2, NULL) == 0) {
+ bond_count++;
+ nval1 = j;
+ nval2 = k2;
+ if (bond_count > 1)
+ return 0; /* Too many bonds */
+ }
+ }
+ }
}
}
if (bond_count == 1) {
return idx;
}
+int
+MoleculeValidatePiConnectionTable(Molecule *mol)
+{
+ Int pass, i, j, k, m;
+
+ if (mol == NULL || mol->pibonds == NULL)
+ return -1; /* No need to process */
+ if (mol->piconnects != NULL)
+ return 0; /* Already valid */
+
+ /* Allocate index table */
+ NewArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), mol->natoms + 1);
+ memset(mol->piconnects, 0, sizeof(Int) * (mol->natoms + 1));
+
+ /* Pass 1: count connections for each atom */
+ /* Pass 2: store connection info */
+ for (pass = 0; pass < 2; pass++) {
+ Int n[2], *ip[2], c[2];
+ for (i = 0; i < mol->npibonds; i++) {
+ AtomConnect *ac;
+ if (mol->pibonds[i * 4 + 2] >= 0)
+ continue; /* Skip angle or dihedral entries */
+ for (j = 0; j < 2; j++) {
+ n[j] = mol->pibonds[i * 4 + j];
+ if (n[j] >= ATOMS_MAX_NUMBER) {
+ ac = &(mol->piatoms[n[j] - ATOMS_MAX_NUMBER].connect);
+ ip[j] = AtomConnectData(ac);
+ c[j] = ac->count;
+ } else if (n[j] >= 0 && n[j] < mol->natoms) {
+ ip[j] = &n[j];
+ c[j] = 1;
+ } else break;
+ }
+ if (j < 2)
+ continue; /* Ignore the invalid entry */
+ for (j = 0; j < c[0]; j++) {
+ for (k = 0; k < c[1]; k++) {
+ Int a1 = ip[0][j];
+ Int a2 = ip[1][k];
+ if (pass == 0) {
+ /* Count */
+ mol->piconnects[a1]++;
+ mol->piconnects[a2]++;
+ } else {
+ /* Store the entry (at the first empty slot) */
+ for (m = mol->piconnects[a1]; m < mol->piconnects[a1 + 1]; m++) {
+ if (mol->piconnects[m] == -1) {
+ mol->piconnects[m] = a2;
+ break;
+ }
+ }
+ for (m = mol->piconnects[a2]; m < mol->piconnects[a2 + 1]; m++) {
+ if (mol->piconnects[m] == -1) {
+ mol->piconnects[m] = a1;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ if (pass == 0) {
+ /* Expand the table, and store the position numbers */
+ m = mol->natoms + 1;
+ for (i = 0; i <= mol->natoms; i++) {
+ j = mol->piconnects[i];
+ mol->piconnects[i] = m;
+ m += j;
+ }
+ AssignArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), m - 1, NULL);
+ for (j = mol->natoms + 1; j < m; j++)
+ mol->piconnects[j] = -1;
+ }
+ }
+ return (mol->npiconnects - mol->natoms - 1); /* Returns the number of entries */
+}
+
+void
+MoleculeInvalidatePiConnectionTable(Molecule *mol)
+{
+ if (mol == NULL)
+ return;
+ if (mol->piconnects != NULL) {
+ free(mol->piconnects);
+ mol->piconnects = NULL;
+ }
+ mol->npiconnects = 0;
+}
+
#pragma mark ====== MO calculation ======
/* Calculate an MO value for a single point. */
where n# is atom index if it is <ATOMS_MAX_NUMBER and
is piatom index + ATOMS_MAX_NUMBER otherwise.
The size of array is 4*npibonds. */
-
+ Int npiconnects; /* Connection table for pi-metal bonds */
+ Int *piconnects; /* The first (natoms + 1) entries are for lookup table, and
+ the connection data (indices of connected atoms) follow.
+ The connected atoms (only by pi-bonds) for atom i are listed in
+ piconnects[piconnects[i]]...piconnects[piconnects[i+1]] */
+
IntGroup *selection;
Int nframes; /* The number of frames (>= 1). This is a cached value, and should be
recalculated from the atoms if it is -1 */
int MoleculeFlushFrames(Molecule *mp);
int MoleculeCalculatePiAtomPosition(Molecule *mol, int idx, Vector *vp);
-
+int MoleculeValidatePiConnectionTable(Molecule *mol);
+void MoleculeInvalidatePiConnectionTable(Molecule *mol);
+
int MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, const Vector *dyp, const Vector *dzp, Int nx, Int ny, Int nz, int (*callback)(double progress, void *ref), void *ref);
int MoleculeGetDefaultMOGrid(Molecule *mp, Int npoints, Vector *op, Vector *xp, Vector *yp, Vector *zp, Int *nx, Int *ny, Int *nz);
const Cube *MoleculeGetCubeAtIndex(Molecule *mp, Int index);