#include <stddef.h>
#include <ctype.h>
#include <math.h>
+#include <float.h>
#include "Missing.h"
#include "Dcd.h"
NewArray(&(dst->connect.u.ptr), &(dst->connect.count), sizeof(Int), src->connect.count);
memmove(dst->connect.u.ptr, src->connect.u.ptr, sizeof(Int) * src->connect.count);
}
+ if (src->anchor != NULL) {
+ dst->anchor = (PiAnchor *)malloc(sizeof(PiAnchor));
+ if (dst->anchor != NULL)
+ memmove(dst->anchor, src->anchor, sizeof(PiAnchor));
+ if (dst->anchor->connect.count > ATOM_CONNECT_LIMIT) {
+ dst->anchor->connect.u.ptr = NULL;
+ dst->anchor->connect.count = 0;
+ NewArray(&(dst->anchor->connect.u.ptr), &(dst->anchor->connect.count), sizeof(Int), src->anchor->connect.count);
+ memmove(dst->anchor->connect.u.ptr, src->anchor->connect.u.ptr, sizeof(Int) * src->anchor->connect.count);
+ }
+ if (dst->anchor->ncoeffs > 0) {
+ NewArray(&(dst->anchor->coeffs), &(dst->anchor->ncoeffs), sizeof(Double), src->anchor->ncoeffs);
+ memmove(dst->anchor->coeffs, src->anchor->coeffs, sizeof(Double) * src->anchor->ncoeffs);
+ }
+ }
return dst;
}
free(bset->moenergies);
if (bset->scfdensities != NULL)
free(bset->scfdensities);
- if (bset->pos != NULL)
- free(bset->pos);
+/* if (bset->pos != NULL)
+ free(bset->pos); */
if (bset->nuccharges != NULL)
free(bset->nuccharges);
if (bset->cubes != NULL) {
return 0;
}
-#if PIATOM
-void
-PiAtomDuplicate(PiAtom *pa, const PiAtom *cpa)
-{
- memmove(pa, cpa, sizeof(PiAtom));
- pa->connect.count = 0;
- AtomConnectResize(&(pa->connect), cpa->connect.count);
- memmove(AtomConnectData(&(pa->connect)), AtomConnectData((AtomConnect *)&(cpa->connect)), sizeof(Int) * cpa->connect.count);
- pa->ncoeffs = 0;
- pa->coeffs = NULL;
- if (cpa->ncoeffs > 0) {
- NewArray(&(pa->coeffs), &(pa->ncoeffs), sizeof(Double), cpa->ncoeffs);
- memmove(pa->coeffs, cpa->coeffs, sizeof(Double) * cpa->ncoeffs);
- }
-}
-
-void
-PiAtomClean(PiAtom *pa)
-{
- AtomConnectResize(&(pa->connect), 0);
- pa->ncoeffs = 0;
- if (pa->coeffs != NULL) {
- free(pa->coeffs);
- pa->coeffs = NULL;
- }
-}
-#endif
-
#pragma mark ====== Accessor types ======
MolEnumerable *
Molecule *
MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
{
+ int i, n;
MoleculeFlushFrames(mp);
MoleculeInitWithAtoms(mp2, mp->atoms, mp->natoms);
if (mp->nbonds > 0) {
NewArray(&(mp2->syms), &(mp2->nsyms), sizeof(Transform), mp->nsyms);
memmove(mp2->syms, mp->syms, sizeof(Transform) * mp2->nsyms);
}
-#if PIATOM
- if (mp->npiatoms > 0) {
- NewArray(&(mp2->piatoms), &(mp2->npiatoms), sizeof(PiAtom), mp->npiatoms);
- for (i = 0; i < mp->npiatoms; i++)
- PiAtomDuplicate(mp2->piatoms + i, mp->piatoms + i);
- }
- if (mp->npibonds > 0) {
- 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);
- }
-#endif
-
-/* mp2->useFlexibleCell = mp->useFlexibleCell; */
+
+ /* mp2->useFlexibleCell = mp->useFlexibleCell; */
if (mp->nframe_cells > 0) {
if (NewArray(&mp2->frame_cells, &mp2->nframe_cells, sizeof(Vector) * 4, mp->nframe_cells) == NULL)
goto error;
memmove(mp2->frame_cells, mp->frame_cells, sizeof(Vector) * 4 * mp->nframe_cells);
}
+ if (mp->nmolprops > 0) {
+ if (NewArray(&mp2->molprops, &mp2->nmolprops, sizeof(MolProp), mp->nmolprops) == NULL)
+ goto error;
+ n = MoleculeGetNumberOfFrames(mp);
+ for (i = 0; i < mp2->nmolprops; i++) {
+ mp2->molprops[i].propname = strdup(mp->molprops[i].propname);
+ mp2->molprops[i].propvals = (Double *)malloc(sizeof(Double) * n);
+ memcpy(mp2->molprops[i].propvals, mp->molprops[i].propvals, sizeof(Double) * n);
+ }
+ }
+
/* FIXME: should bset (basis set info) and elpot be duplicated or not? */
if (mp->par != NULL)
void
MoleculeClear(Molecule *mp)
{
+ int i;
if (mp == NULL)
return;
if (mp->arena != NULL) {
ParameterRelease(mp->par);
mp->par = NULL;
}
- if (mp->bset != NULL) {
- BasisSetRelease(mp->bset);
- mp->bset = NULL;
- }
if (mp->atoms != NULL) {
- int i;
for (i = 0; i < mp->natoms; i++)
AtomClean(mp->atoms + i);
free(mp->atoms);
mp->syms = NULL;
mp->nsyms = 0;
}
-#if PIATOM
- if (mp->piatoms != NULL) {
- for (i = 0; i < mp->npiatoms; i++) {
- PiAtomClean(mp->piatoms + i);
- }
- free(mp->piatoms);
- mp->piatoms = NULL;
- mp->npiatoms = 0;
- }
- if (mp->pibonds != NULL) {
- free(mp->pibonds);
- mp->pibonds = NULL;
- mp->npibonds = 0;
- }
- if (mp->piconnects != NULL) {
- free(mp->piconnects);
- mp->piconnects = NULL;
- mp->npiconnects = 0;
- }
-#endif
if (mp->selection != NULL) {
IntGroupRelease(mp->selection);
mp->selection = NULL;
BasisSetRelease(mp->bset);
mp->bset = NULL;
}
+ if (mp->mcube != NULL) {
+ MoleculeDeallocateMCube(mp->mcube);
+ mp->mcube = NULL;
+ }
+ if (mp->molprops != NULL) {
+ for (i = 0; i < mp->nmolprops; i++) {
+ free(mp->molprops[i].propname);
+ free(mp->molprops[i].propvals);
+ }
+ free(mp->molprops);
+ mp->molprops = NULL;
+ mp->nmolprops = 0;
+ }
if (mp->par != NULL) {
ParameterRelease(mp->par);
mp->par = NULL;
int ibuf[12];
Int iibuf[4];
double dbuf[12];
- int mview_ibuf[16];
- float mview_fbuf[8];
+ int mview_ibuf[18];
+ double mview_dbuf[10];
char cbuf[12][8];
const char **pp;
char *bufp, *valp, *comp;
s_append_asprintf(errbuf, "Cannot open file");
return 1;
}
- for (i = 0; i < 8; i++)
- mview_fbuf[i] = kUndefined;
- for (i = 0; i < 16; i++)
+ for (i = 0; i < 10; i++)
+ mview_dbuf[i] = kUndefined;
+ for (i = 0; i < 18; i++)
mview_ibuf[i] = kUndefined;
/* flockfile(fp); */
lineNumber = 0;
i++;
}
continue;
+ } else if (strcmp(buf, "!:uff_types") == 0) {
+ i = 0;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ /* idx uff_type */
+ if (sscanf(buf, "%d %6s", &ibuf[0], cbuf[0]) < 2) {
+ s_append_asprintf(errbuf, "line %d: uff type info cannot be read for atom %d", lineNumber, i + 1);
+ goto err_exit;
+ }
+ if (i >= mp->natoms) {
+ s_append_asprintf(errbuf, "line %d: too many uff type info\n", lineNumber);
+ goto err_exit;
+ }
+ ap = ATOM_AT_INDEX(mp->atoms, i);
+ strncpy(ap->uff_type, cbuf[0], 5);
+ ap->uff_type[5] = 0;
+ i++;
+ }
} else if (strcmp(buf, "!:mm_exclude") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
i++;
}
continue;
+ } else if (strcmp(buf, "!:pi_anchor") == 0) {
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ /* idx count */
+ if ((j = sscanf(buf, "%d %d", &ibuf[0], &ibuf[1])) < 2) {
+ s_append_asprintf(errbuf, "line %d: bad format for pi_anchor", lineNumber);
+ goto err_exit;
+ }
+ i = ibuf[0];
+ ap = ATOM_AT_INDEX(mp->atoms, i);
+ if (ap->anchor != NULL) {
+ s_append_asprintf(errbuf, "line %d: warning: duplicate pi_anchor entry", lineNumber);
+ AtomConnectResize(&ap->anchor->connect, 0);
+ free(ap->anchor->coeffs);
+ free(ap->anchor);
+ }
+ ap->anchor = (PiAnchor *)calloc(sizeof(PiAnchor), 1);
+ if (ibuf[1] < 2 || ibuf[1] >= mp->natoms) {
+ s_append_asprintf(errbuf, "line %d: bad number of components for pi_anchor", lineNumber);
+ goto err_exit;
+ }
+ AtomConnectResize(&ap->anchor->connect, ibuf[1]);
+ ip = AtomConnectData(&ap->anchor->connect);
+ NewArray(&ap->anchor->coeffs, &ap->anchor->ncoeffs, sizeof(Double), ibuf[1]);
+ j = ibuf[1];
+ for (i = 0; i < j; i++) {
+ if (ReadLine(buf, sizeof buf, fp, &lineNumber) <= 0) {
+ s_append_asprintf(errbuf, "line %d: unexpected end of file while reading pi_anchors", lineNumber);
+ goto err_exit;
+ }
+ if (sscanf(buf, "%d %lf", &ibuf[0], &dbuf[0]) < 2) {
+ s_append_asprintf(errbuf, "line %d: bad format for pi_anchor", lineNumber);
+ goto err_exit;
+ }
+ if (ibuf[0] < 0 || ibuf[0] >= mp->natoms) {
+ s_append_asprintf(errbuf, "line %d: atom index out of range", lineNumber);
+ goto err_exit;
+ }
+ if (dbuf[0] <= 0.0) {
+ s_append_asprintf(errbuf, "line %d: the pi anchor weights should be positive", lineNumber);
+ goto err_exit;
+ }
+ ip[i] = ibuf[0];
+ ap->anchor->coeffs[i] = dbuf[0];
+ }
+ }
+ continue;
} else if (strcmp(buf, "!:positions") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
}
}
continue;
+ } else if (strcmp(buf, "!:bond_orders") == 0) {
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ /* b1 b2 b3 b4 */
+ i = sscanf(buf, "%lf %lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2], &dbuf[3]);
+ if (i == 0) {
+ s_append_asprintf(errbuf, "line %d: bad bond order format", lineNumber);
+ goto err_exit;
+ }
+ for (j = 0; j < i; j++) {
+ AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbondOrders, &dbuf[j]);
+ }
+ }
+ if (mp->nbondOrders > mp->nbonds) {
+ s_append_asprintf(errbuf, "line %d: warning: the number of bond order info (%d) exceeds number of bonds (%d) - ignoring excess info\n", lineNumber, mp->nbondOrders, mp->nbonds);
+ nwarnings++;
+ mp->nbondOrders = mp->nbonds;
+ } else if (mp->nbondOrders < mp->nbonds) {
+ s_append_asprintf(errbuf, "line %d: warning: the number of bond order info (%d) is less than number of bonds (%d)\n", lineNumber, mp->nbondOrders, mp->nbonds);
+ nwarnings++;
+ j = mp->nbondOrders;
+ AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbonds - 1, NULL);
+ for (i = j; i < mp->nbonds; i++)
+ mp->bondOrders[i] = 0.0;
+ }
+ continue;
+
} else if (strcmp(buf, "!:angles") == 0) {
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
if (buf[0] == '!')
if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2]) {
s_append_asprintf(errbuf, "line %d: warning: bad angle specification (%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2]);
nwarnings++;
- } else if (MoleculeLookupBond(mp, iibuf[0], iibuf[1]) < 0 || MoleculeLookupBond(mp, iibuf[1], iibuf[2]) < 0) {
+ } else if (MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0) {
s_append_asprintf(errbuf, "line %d: warning: angle with non-bonded atoms (%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2]);
nwarnings++;
} else if (MoleculeLookupAngle(mp, iibuf[0], iibuf[1], iibuf[2]) >= 0) {
if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[3] < 0 || iibuf[3] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2] || iibuf[2] == iibuf[3] || iibuf[0] == iibuf[2] || iibuf[1] == iibuf[3] || iibuf[0] == iibuf[3]) {
s_append_asprintf(errbuf, "line %d: warning: bad dihedral specification (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
nwarnings++;
- } else if (MoleculeLookupBond(mp, iibuf[0], iibuf[1]) < 0 || MoleculeLookupBond(mp, iibuf[1], iibuf[2]) < 0 || MoleculeLookupBond(mp, iibuf[2], iibuf[3]) < 0) {
+ } else if (MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
s_append_asprintf(errbuf, "line %d: warning: dihedral with non-bonded atoms (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
nwarnings++;
} else if (MoleculeLookupDihedral(mp, iibuf[0], iibuf[1], iibuf[2], iibuf[3]) >= 0) {
if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[3] < 0 || iibuf[3] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2] || iibuf[2] == iibuf[3] || iibuf[0] == iibuf[2] || iibuf[1] == iibuf[3] || iibuf[0] == iibuf[3]) {
s_append_asprintf(errbuf, "line %d: warning: bad improper specification (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
nwarnings++;
- } else if (MoleculeLookupBond(mp, iibuf[0], iibuf[2]) < 0 || MoleculeLookupBond(mp, iibuf[1], iibuf[2]) < 0 || MoleculeLookupBond(mp, iibuf[2], iibuf[3]) < 0) {
+ } else if (MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[1]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
s_append_asprintf(errbuf, "line %d: warning: improper with non-bonded atoms (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
nwarnings++;
} else if (MoleculeLookupImproper(mp, iibuf[0], iibuf[1], iibuf[2], iibuf[3]) >= 0) {
}
}
continue;
-#if PIATOM
- } else if (strcmp(buf, "!:pi_atoms") == 0) {
- PiAtom *pp;
- while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
- if (buf[0] == '!')
- continue;
- if (buf[0] == '\n')
- break;
- if (sscanf(buf, "%6s %6s %d", cbuf[0], cbuf[1], &ibuf[0]) < 3) {
- s_append_asprintf(errbuf, "line %d: pi atoms info cannot be read", lineNumber);
- goto err_exit;
- }
- pp = (PiAtom *)AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, NULL);
- memset(pp, 0, sizeof(PiAtom));
- strncpy(pp->aname, cbuf[0], 4);
- pp->type = AtomTypeEncodeToUInt(cbuf[1]);
- if (ibuf[0] <= 0)
- continue;
- AtomConnectResize(&pp->connect, ibuf[0]);
- ip = AtomConnectData(&pp->connect);
- NewArray(&pp->coeffs, &pp->ncoeffs, sizeof(Double), ibuf[0]);
- for (i = 0; i < ibuf[0]; i++) {
- if (ReadLine(buf, sizeof buf, fp, &lineNumber) <= 0) {
- s_append_asprintf(errbuf, "line %d: unexpected end of file during reading pi atoms info", lineNumber);
- goto err_exit;
- }
- if (sscanf(buf, "%d %lf", &ibuf[1], &dbuf[0]) < 2) {
- s_append_asprintf(errbuf, "line %d: bad format during pi atoms info", lineNumber);
- goto err_exit;
- }
- ip[i] = ibuf[1];
- pp->coeffs[i] = dbuf[0];
- }
- }
- continue;
- } else if (strcmp(buf, "!:pi_atom_constructs") == 0) {
- while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
- if (buf[0] == '!')
- continue;
- if (buf[0] == '\n')
- break;
- /* a1 b1 c1 d1 a2 b2 c2 d2 */
- i = sscanf(buf, "%d %d %d %d %d %d %d %d", &ibuf[0], &ibuf[1], &ibuf[2], &ibuf[3], &ibuf[4], &ibuf[5], &ibuf[6], &ibuf[7]);
- if (i == 0 || i % 4 != 0)
- goto pi_atom_constructs_bad_format;
- for (j = 0; j < i; j++) {
- if (ibuf[j] <= -2) {
- ibuf[j] = (-ibuf[j] - 2) + ATOMS_MAX_NUMBER;
- if (ibuf[j] - ATOMS_MAX_NUMBER >= mp->npiatoms)
- goto pi_atom_constructs_bad_format;
- } else if (ibuf[j] >= mp->natoms) {
- goto pi_atom_constructs_bad_format;
- }
- if (j % 4 == 3) {
- AssignArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, mp->npibonds, &ibuf[j - 3]);
- }
- }
- }
- continue;
- pi_atom_constructs_bad_format:
- s_append_asprintf(errbuf, "line %d: bad format in pi_atom_constructs", lineNumber);
- goto err_exit;
-#endif /* PIATOM */
} else if (strcmp(buf, "!:anisotropic_thermal_parameters") == 0) {
i = 0;
while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
|| (strcmp(comp, "cutoff") == 0 && (dp = &arena->cutoff) != NULL)
|| (strcmp(comp, "electro_cutoff") == 0 && (dp = &arena->electro_cutoff) != NULL)
|| (strcmp(comp, "pairlist_distance") == 0 && (dp = &arena->pairlist_distance) != NULL)
+ || (strcmp(comp, "switch_distance") == 0 && (dp = &arena->switch_distance) != NULL)
|| (strcmp(comp, "temperature") == 0 && (dp = &arena->temperature) != NULL)
|| (strcmp(comp, "andersen_coupling") == 0 && (dp = &arena->andersen_thermo_coupling) != NULL)
|| (strcmp(comp, "dielectric") == 0 && (dp = &arena->dielectric) != NULL)
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;
continue;
if (buf[0] == '\n')
break;
+ if (mp->mview == NULL || mp->mview->track == NULL)
+ continue; /* Skip (this should not happen though) */
/* scale; trx try trz; theta_deg x y z */
- if ((i == 0 && sscanf(buf, "%f", &mview_fbuf[0]) < 1)
- || (i == 1 && sscanf(buf, "%f %f %f",
- &mview_fbuf[1], &mview_fbuf[2], &mview_fbuf[3]) < 3)
- || (i == 2 && sscanf(buf, "%f %f %f %f",
- &mview_fbuf[4], &mview_fbuf[5], &mview_fbuf[6], &mview_fbuf[7]) < 4)) {
+ if ((i == 0 && sscanf(buf, "%lf", &dbuf[0]) < 1)
+ || (i == 1 && sscanf(buf, "%lf %lf %lf",
+ &dbuf[1], &dbuf[2], &dbuf[3]) < 3)
+ || (i == 2 && sscanf(buf, "%lf %lf %lf %lf",
+ &dbuf[4], &dbuf[5], &dbuf[6], &dbuf[7]) < 4)) {
s_append_asprintf(errbuf, "line %d: bad trackball format", lineNumber);
goto err_exit;
}
+ if (i == 0)
+ TrackballSetScale(mp->mview->track, dbuf[0]);
+ else if (i == 1)
+ TrackballSetTranslate(mp->mview->track, dbuf + 1);
+ else if (i == 2)
+ TrackballSetRotate(mp->mview->track, dbuf + 4);
i++;
}
continue;
continue;
if (buf[0] == '\n')
break;
+ if (mp->mview == NULL)
+ continue; /* Skip (this should not happen, though) */
bufp = buf;
comp = strsep(&bufp, " \t");
if (bufp != NULL) {
bufp++;
valp = strsep(&bufp, "\n");
} else valp = NULL;
- /* In the following, the redundant "!= NULL" is to suppress suprious warning */
- if ((strcmp(comp, "show_unit_cell") == 0 && (i = 1))
- || (strcmp(comp, "show_periodic_box") == 0 && (i = 2))
- || (strcmp(comp, "show_expanded_atoms") == 0 && (i = 3))
- || (strcmp(comp, "show_ellipsoids") == 0 && (i = 4))
- || (strcmp(comp, "show_hydrogens") == 0 && (i = 5))
- || (strcmp(comp, "show_dummy_atoms") == 0 && (i = 6))
- || (strcmp(comp, "show_rotation_center") == 0 && (i = 7))
- || (strcmp(comp, "show_graphite_flag") == 0 && (i = 8))
- || (strcmp(comp, "show_periodic_image_flag") == 0 && (i = 9))
- || (strcmp(comp, "show_graphite") == 0 && (i = 10))) {
- mview_ibuf[i - 1] = atoi(valp);
- } else if (strcmp(comp, "show_periodic_image") == 0) {
- sscanf(valp, "%d %d %d %d %d %d",
- &mview_ibuf[10], &mview_ibuf[11], &mview_ibuf[12],
- &mview_ibuf[13], &mview_ibuf[14], &mview_ibuf[15]);
+ if (strcmp(comp, "show_unit_cell") == 0)
+ mp->mview->showUnitCell = atoi(valp);
+ else if (strcmp(comp, "show_periodic_box") == 0)
+ mp->mview->showPeriodicBox = atoi(valp);
+ else if (strcmp(comp, "show_expanded_atoms") == 0)
+ mp->mview->showExpandedAtoms = atoi(valp);
+ else if (strcmp(comp, "show_ellipsoids") == 0)
+ mp->mview->showEllipsoids = atoi(valp);
+ else if (strcmp(comp, "show_hydrogens") == 0)
+ mp->mview->showHydrogens = atoi(valp);
+ else if (strcmp(comp, "show_dummy_atoms") == 0)
+ mp->mview->showDummyAtoms = atoi(valp);
+ else if (strcmp(comp, "show_rotation_center") == 0)
+ mp->mview->showRotationCenter = atoi(valp);
+ else if (strcmp(comp, "show_graphite_flag") == 0)
+ mp->mview->showGraphiteFlag = atoi(valp);
+ else if (strcmp(comp, "show_periodic_image_flag") == 0)
+ mp->mview->showPeriodicImageFlag = atoi(valp);
+ else if (strcmp(comp, "show_graphite") == 0)
+ mp->mview->showGraphite = atoi(valp);
+ else if (strcmp(comp, "show_expanded_atoms") == 0)
+ mp->mview->showExpandedAtoms = atoi(valp);
+ else if (strcmp(comp, "atom_resolution") == 0 && (i = atoi(valp)) >= 6)
+ mp->mview->atomResolution = i;
+ else if (strcmp(comp, "bond_resolution") == 0 && (i = atoi(valp)) >= 4)
+ mp->mview->bondResolution = i;
+ else if (strcmp(comp, "atom_radius") == 0)
+ mp->mview->atomRadius = strtod(valp, NULL);
+ else if (strcmp(comp, "bond_radius") == 0)
+ mp->mview->bondRadius = strtod(valp, NULL);
+ else if (strcmp(comp, "show_periodic_image") == 0) {
+ sscanf(valp, "%d %d %d %d %d %d", &ibuf[0], &ibuf[1], &ibuf[2], &ibuf[3], &ibuf[4], &ibuf[5]);
+ for (i = 0; i < 6; i++)
+ mp->mview->showPeriodicImage[i] = ibuf[i];
+ }
+ }
+ continue;
+ } else if (strcmp(buf, "!:property") == 0) {
+ char dec[1024];
+ i = 0;
+ bufp = buf + 13;
+ while (*bufp != 0 && *bufp != '\n' && bufp < (buf + sizeof buf - 3)) {
+ if (*bufp == '%') {
+ dec[i] = bufp[1];
+ dec[i + 1] = bufp[2];
+ dec[i + 2] = 0;
+ dec[i++] = strtol(dec, NULL, 16);
+ bufp += 3;
+ } else {
+ dec[i++] = *bufp++;
+ }
+ if (i >= 1000)
+ break;
+ }
+ if (i == 0)
+ continue;
+ dec[i] = 0;
+ i = MoleculeCreateProperty(mp, dec);
+ if (i < 0) {
+ s_append_asprintf(errbuf, "line %d: warning: duplicate molecular property %s - ignored\n", lineNumber, dec);
+ nwarnings++;
+ continue;
+ }
+ j = 0;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (j >= nframes) {
+ s_append_asprintf(errbuf, "line %d: warning: too many molecular property %s - ignored\n", lineNumber, dec);
+ nwarnings++;
+ break;
+ }
+ dbuf[0] = strtod(buf, NULL);
+ mp->molprops[i].propvals[j] = dbuf[0];
+ j++;
+ }
+ continue;
+ } else if (strcmp(buf, "!:gaussian_primitives") == 0) {
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ /* sym nprims a_idx */
+ if (sscanf(buf, "%6s %d %d", cbuf[0], &ibuf[0], &ibuf[1]) < 3) {
+ s_append_asprintf(errbuf, "line %d: the gaussian primitive info cannot be read", lineNumber);
+ goto err_exit;
+ }
+ if (strcasecmp(cbuf[0], "S") == 0) {
+ ibuf[2] = 0;
+ } else if (strcasecmp(cbuf[0], "P") == 0) {
+ ibuf[2] = 1;
+ } else if (strcasecmp(cbuf[0], "SP") == 0) {
+ ibuf[2] = -1;
+ } else if (strcasecmp(cbuf[0], "D") == 0) {
+ ibuf[2] = 2;
+ } else if (strcasecmp(cbuf[0], "D5") == 0) {
+ ibuf[2] = -2;
+ } else if (strcasecmp(cbuf[0], "F") == 0) {
+ ibuf[2] = 3;
+ } else if (strcasecmp(cbuf[0], "F7") == 0) {
+ ibuf[2] = -3;
+ } else if (strcasecmp(cbuf[0], "G") == 0) {
+ ibuf[2] = 4;
+ } else if (strcasecmp(cbuf[0], "G9") == 0) {
+ ibuf[2] = -4;
+ } else {
+ s_append_asprintf(errbuf, "line %d: the gaussian primitive type %s is unknown", lineNumber, cbuf[0]);
+ goto err_exit;
+ }
+ if (ibuf[0] <= 0) {
+ s_append_asprintf(errbuf, "line %d: the number of primitive (%d) must be positive", lineNumber, ibuf[0]);
+ goto err_exit;
+ }
+ if (ibuf[1] < 0 || ibuf[1] >= mp->natoms) {
+ s_append_asprintf(errbuf, "line %d: the atom index (%d) is out of range", lineNumber, ibuf[1]);
+ goto err_exit;
+ }
+ MoleculeAddGaussianOrbitalShell(mp, ibuf[1], ibuf[2], ibuf[0], 0);
+ i = ibuf[0];
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (sscanf(buf, "%lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2]) < 3) {
+ s_append_asprintf(errbuf, "line %d: cannot read gaussian primitive coefficients", lineNumber);
+ goto err_exit;
+ }
+ MoleculeAddGaussianPrimitiveCoefficients(mp, dbuf[0], dbuf[1], dbuf[2]);
+ if (--i == 0)
+ break;
+ }
+ if (buf[0] == '\n')
+ break;
+ }
+ continue;
+ } else if (strcmp(buf, "!:mo_info") == 0) {
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (sscanf(buf, "%6s %d %d", cbuf[0], &ibuf[0], &ibuf[1]) < 3) {
+ s_append_asprintf(errbuf, "line %d: the MO info cannot be correctly read", lineNumber);
+ goto err_exit;
+ }
+ if (strcasecmp(cbuf[0], "RHF") == 0) {
+ ibuf[2] = 1;
+ } else if (strcasecmp(cbuf[0], "ROHF") == 0) {
+ ibuf[2] = 2;
+ } else if (strcasecmp(cbuf[0], "UHF") == 0) {
+ ibuf[2] = 0;
+ } else {
+ s_append_asprintf(errbuf, "line %d: unknown HF type: %s", lineNumber, cbuf[0]);
+ goto err_exit;
+ }
+ if (ibuf[0] < 0 || ibuf[1] < 0) {
+ s_append_asprintf(errbuf, "line %d: incorrect number of electrons", lineNumber);
+ goto err_exit;
+ }
+ MoleculeSetMOInfo(mp, ibuf[2], ibuf[0], ibuf[1]);
+ }
+ continue;
+ } else if (strcmp(buf, "!:mo_coefficients") == 0) {
+ if (mp->bset == NULL || mp->bset->nshells == 0) {
+ s_append_asprintf(errbuf, "line %d: the :gaussian_primitive section must come before :mo_coefficients", lineNumber);
+ goto err_exit;
+ }
+ /* Count the number of components */
+ dp = (Double *)malloc(sizeof(Double) * mp->bset->ncomps);
+ i = 1;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (sscanf(buf, "MO %d %lf", &ibuf[0], &dbuf[6]) < 2) {
+ s_append_asprintf(errbuf, "line %d: cannot read the MO index or energy", lineNumber);
+ goto err_exit;
+ }
+ if (ibuf[0] != i) {
+ s_append_asprintf(errbuf, "line %d: the MO index (%d) must be in ascending order", lineNumber, ibuf[0]);
+ goto err_exit;
+ }
+ i = 0;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ j = sscanf(buf, "%lf %lf %lf %lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2], &dbuf[3], &dbuf[4], &dbuf[5]);
+ if (j == 0) {
+ s_append_asprintf(errbuf, "line %d: cannot read the MO coefficients", lineNumber);
+ goto err_exit;
+ }
+ for (k = 0; k < j; k++, i++) {
+ if (i >= mp->bset->ncomps) {
+ s_append_asprintf(errbuf, "line %d: too many MO coefficients", lineNumber);
+ goto err_exit;
+ }
+ dp[i] = dbuf[k];
+ }
+ if (i >= mp->bset->ncomps)
+ break;
}
+ i = MoleculeSetMOCoefficients(mp, ibuf[0], dbuf[6], mp->bset->ncomps, dp);
+ if (i != 0) {
+ s_append_asprintf(errbuf, "line %d: cannot set MO coefficients", lineNumber);
+ goto err_exit;
+ }
+ i = ibuf[0] + 1; /* For next entry */
+ }
+ continue;
+ } else if (strcmp(buf, "!:graphics") == 0) {
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ MainViewGraphic *gp = NULL;
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (mp->mview == NULL)
+ continue; /* Skip */
+ redo:
+ if (strcmp(buf, "line\n") == 0) {
+ ibuf[0] = kMainViewGraphicLine;
+ } else if (strcmp(buf, "poly\n") == 0) {
+ ibuf[0] = kMainViewGraphicPoly;
+ } else if (strcmp(buf, "cylinder\n") == 0) {
+ ibuf[0] = kMainViewGraphicCylinder;
+ } else if (strcmp(buf, "cone\n") == 0) {
+ ibuf[0] = kMainViewGraphicCone;
+ } else if (strcmp(buf, "ellipsoid\n") == 0) {
+ ibuf[0] = kMainViewGraphicEllipsoid;
+ } else {
+ continue; /* Skip */
+ }
+ gp = (MainViewGraphic *)calloc(sizeof(MainViewGraphic), 1);
+ gp->kind = ibuf[0];
+ i = 0;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ if (buf[0] == '!')
+ continue;
+ if (buf[0] == '\n')
+ break;
+ if (i == 0) {
+ if (sscanf(buf, "%d %d", &ibuf[0], &ibuf[1]) < 2) {
+ s_append_asprintf(errbuf, "line %d: the closed/visible flags cannot be read for graphic object", lineNumber);
+ goto err_exit;
+ }
+ gp->closed = ibuf[0];
+ gp->visible = ibuf[1];
+ } else if (i == 1) {
+ if (sscanf(buf, "%lf %lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2], &dbuf[3]) < 4) {
+ s_append_asprintf(errbuf, "line %d: the color cannot be read for graphic object", lineNumber);
+ goto err_exit;
+ }
+ for (j = 0; j < 4; j++)
+ gp->rgba[j] = dbuf[j];
+ } else if (i == 2) {
+ j = atoi(buf);
+ if (j < 0) {
+ s_append_asprintf(errbuf, "line %d: the number of control points must be non-negative", lineNumber);
+ goto err_exit;
+ }
+ if (j > 0)
+ NewArray(&gp->points, &gp->npoints, sizeof(GLfloat) * 3, j);
+ } else if (i >= 3 && i < gp->npoints + 3) {
+ if (sscanf(buf, "%lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2]) < 3) {
+ s_append_asprintf(errbuf, "line %d: the control point cannot be read for graphic object", lineNumber);
+ goto err_exit;
+ }
+ j = (i - 3) * 3;
+ gp->points[j++] = dbuf[0];
+ gp->points[j++] = dbuf[1];
+ gp->points[j] = dbuf[2];
+ } else if (i == gp->npoints + 3) {
+ j = atoi(buf);
+ if (j < 0) {
+ s_append_asprintf(errbuf, "line %d: the number of normals must be non-negative", lineNumber);
+ goto err_exit;
+ }
+ if (j > 0)
+ NewArray(&gp->normals, &gp->nnormals, sizeof(GLfloat) * 3, j);
+ } else if (i >= gp->npoints + 4 && i < gp->npoints + gp->nnormals + 4) {
+ if (sscanf(buf, "%lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2]) < 3) {
+ s_append_asprintf(errbuf, "line %d: the normal vector cannot be read for graphic object", lineNumber);
+ goto err_exit;
+ }
+ j = (i - gp->npoints - 4) * 3;
+ gp->normals[j++] = dbuf[0];
+ gp->normals[j++] = dbuf[1];
+ gp->normals[j] = dbuf[2];
+ } else break;
+ i++;
+ }
+ MainView_insertGraphic(mp->mview, -1, gp);
+ free(gp);
+ if (buf[0] == '\n' || buf[0] == 0)
+ break;
+ goto redo;
+ }
+ continue;
+ } else if (strncmp(buf, "!:@", 3) == 0) {
+ /* Plug-in implemented in the ruby world */
+ Int stringLen;
+ char *stringBuf, *returnString;
+ i = strlen(buf);
+ NewArray(&stringBuf, &stringLen, sizeof(char), i + 1);
+ strcpy(stringBuf, buf);
+ k = lineNumber;
+ while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+ /* The comment lines are _not_ skipped */
+ if (buf[0] == '\n')
+ break;
+ j = strlen(buf);
+ AssignArray(&stringBuf, &stringLen, sizeof(char), i + j, NULL);
+ strncpy(stringBuf + i, buf, j);
+ i += j;
+ }
+ if (MolActionCreateAndPerform(mp, SCRIPT_ACTION("si;s"),
+ "proc { |i| loadmbsf_plugin(i) rescue \"line #{i}: #{$i.to_s}\" }",
+ stringBuf, k, &returnString) != 0) {
+ s_append_asprintf(errbuf, "line %d: cannot invoke Ruby plugin", lineNumber);
+ goto err_exit;
+ } else if (returnString[0] != 0) {
+ s_append_asprintf(errbuf, "%s", returnString);
+ goto err_exit;
}
+ free(stringBuf);
continue;
}
/* Unknown sections are silently ignored */
md_arena_set_molecule(mp->arena, mp);
fclose(fp);
- if (mp->mview != NULL) {
+
+/* if (mp->mview != NULL) {
if (mview_ibuf[0] != kUndefined)
mp->mview->showUnitCell = mview_ibuf[0];
if (mview_ibuf[1] != kUndefined)
mp->mview->showPeriodicImageFlag = mview_ibuf[8];
if (mview_ibuf[9] != kUndefined)
mp->mview->showGraphite = mview_ibuf[9];
+ if (mview_ibuf[10] != kUndefined && mview_ibuf[10] >= 6)
+ mp->mview->atomResolution = mview_ibuf[10];
+ if (mview_ibuf[11] != kUndefined && mview_ibuf[11] >= 4)
+ mp->mview->bondResolution = mview_ibuf[11];
for (i = 0; i < 6; i++) {
- if (mview_ibuf[10 + i] != kUndefined)
- mp->mview->showPeriodicImage[i] = mview_ibuf[10 + i];
+ if (mview_ibuf[12 + i] != kUndefined)
+ mp->mview->showPeriodicImage[i] = mview_ibuf[12 + i];
}
+ if (mview_dbuf[8] != kUndefined)
+ mp->mview->atomRadius = mview_dbuf[8];
+ if (mview_dbuf[9] != kUndefined)
+ mp->mview->bondRadius = mview_dbuf[9];
if (mp->mview->track != NULL) {
- if (mview_fbuf[0] != kUndefined)
- TrackballSetScale(mp->mview->track, mview_fbuf[0]);
- if (mview_fbuf[1] != kUndefined)
- TrackballSetTranslate(mp->mview->track, mview_fbuf + 1);
- if (mview_fbuf[4] != kUndefined)
- TrackballSetRotate(mp->mview->track, mview_fbuf + 4);
+ if (mview_dbuf[0] != kUndefined)
+ TrackballSetScale(mp->mview->track, mview_dbuf[0]);
+ if (mview_dbuf[1] != kUndefined)
+ TrackballSetTranslate(mp->mview->track, mview_dbuf + 1);
+ if (mview_dbuf[4] != kUndefined)
+ TrackballSetRotate(mp->mview->track, mview_dbuf + 4);
}
}
+*/
return 0;
frames = (Vector *)realloc(frames, size);
if (frames == NULL)
goto panic;
- #if 0
- if (fn == 1) {
- /* Copy the coordinates of the first frame */
- for (i = 0; i < mp->natoms; i++) {
- ap = ATOM_AT_INDEX(mp->atoms, i);
- frames[i] = ap->r;
- }
- }
- /* Copy the coordinates of the last frame to the newly created frame */
- memmove(frames + sizeof(Vector) * mp->natoms * fn, frames + sizeof(Vector) * mp->natoms * (fn - 1), sizeof(Vector) * mp->natoms);
- #endif
}
/* Read coordinates */
for (i = 0; i < mp->natoms; i++) {
if (fn > 1) {
for (i = 0; i < mp->natoms; i++) {
ap = ATOM_AT_INDEX(mp->atoms, i);
- ap->frames = (Vector *)malloc(sizeof(Vector) * fn);
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), fn);
if (ap->frames == NULL)
goto panic;
- ap->nframes = fn;
for (j = 0; j < fn; j++)
ap->frames[j] = frames[mp->natoms * j + i];
}
}
}
fclose(fp);
- MoleculeGuessBonds(mp, 1.2, &nbonds, &bonds);
+ MoleculeGuessBonds(mp, 0.0, &nbonds, &bonds);
if (nbonds > 0) {
MoleculeAddBonds(mp, nbonds, bonds, NULL, 1);
free(bonds);
sMoleculeGenerateSymopWithTransform(mp, tr_inv, 0);
}
- MoleculeGuessBonds(mp, 1.2, &nbonds, &bonds);
+ MoleculeGuessBonds(mp, 0.0, &nbonds, &bonds);
if (nbonds > 0) {
MoleculeAddBonds(mp, nbonds, bonds, NULL, 1);
free(bonds);
/* Add one gaussian orbital shell information (not undoable) */
int
-MoleculeAddGaussianOrbitalShell(Molecule *mol, Int sym, Int nprims, Int a_idx)
+MoleculeAddGaussianOrbitalShell(Molecule *mol, Int a_idx, Int sym, Int nprims, Int add_exp)
{
BasisSet *bset;
ShellInfo *shellp;
case -1: shellp->sym = kGTOType_SP; shellp->ncomp = 4; break;
case 2: shellp->sym = kGTOType_D; shellp->ncomp = 6; break;
case -2: shellp->sym = kGTOType_D5; shellp->ncomp = 5; break;
- /* TODO: Support F/F7 type orbitals */
- /* case 3: sp->sym = kGTOtype_F; sp->ncomp = 10; break;
- case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break; */
+ case 3: shellp->sym = kGTOType_F; shellp->ncomp = 10; break;
+ case -3: shellp->sym = kGTOType_F7; shellp->ncomp = 7; break;
+ case 4: shellp->sym = kGTOType_G; shellp->ncomp = 15; break;
+ case -4: shellp->sym = kGTOType_G9; shellp->ncomp = 9; break;
default:
return -3; /* Unsupported shell type */
}
shellp->m_idx = 0;
shellp->p_idx = 0;
}
+ shellp->add_exp = add_exp;
+ /* Update the number of components (if not yet determined) */
+ if (bset->ncomps < shellp->m_idx + shellp->ncomp)
+ bset->ncomps = shellp->m_idx + shellp->ncomp;
return 0;
}
return 0;
}
-/* Set MO coefficients for idx-th MO */
+/* Get the shell information from the component index */
+/* The outLabel must have space for at least 23 non-Null characters */
+int
+MoleculeGetGaussianComponentInfo(Molecule *mol, Int comp_idx, Int *outAtomIdx, char *outLabel, Int *outShellIdx)
+{
+ BasisSet *bset;
+ ShellInfo *shellp;
+ int si;
+ if (mol == NULL || (bset = mol->bset) == NULL)
+ return -1; /* No basis set info */
+ if (comp_idx < 0 || comp_idx >= bset->ncomps)
+ return -2; /* Component index out of range */
+ for (si = 0, shellp = bset->shells; si < bset->nshells; si++, shellp++) {
+ if (comp_idx >= shellp->ncomp) {
+ comp_idx -= shellp->ncomp;
+ continue;
+ } else {
+ static const char *type_p = "xyz";
+ static const char *type_d = "xxyyzzxyxzyz";
+ static const char *type_d5[] = {"xy","yz","zz", "xz", "xx-yy"};
+ static const char *type_f = "xxxyyyzzzxxyxxzxyyyyzxzzyzzxyz";
+ static const char *type_f7[] = {"x3-3xy2", "x2z-y2z", "x(5z2-r2)", "z(5z2-3r2)", "y(5z2-r2)", "xyz", "3x2y-y3"};
+ static const char *type_g[] = {"x4", "y4", "z4", "x3y", "x3z", "xy3", "y3z", "xz3", "yz3", "x2y2", "x2z2", "y2z2", "x2yz", "x2yz", "xyz2"};
+ static const char *type_g9[] = {"x4+y4-6x2y2", "xz(x2-3y2)", "(x2-y2)(7z2-r2)", "xz(7z2-3r2)", "35z4-30z2r2+3r4", "yz(7z2-3r2)", "xy(7z2-r2)", "yz(3x2-y2)", "xy(x2-y2)"};
+ *outAtomIdx = shellp->a_idx;
+ *outShellIdx = si;
+ switch (shellp->sym) {
+ case kGTOType_S:
+ strcpy(outLabel, "S");
+ break;
+ case kGTOType_P:
+ outLabel[0] = 'P';
+ outLabel[1] = type_p[comp_idx];
+ outLabel[2] = 0;
+ break;
+ case kGTOType_SP:
+ if (comp_idx == 0)
+ strcpy(outLabel, "S");
+ else {
+ outLabel[0] = 'P';
+ outLabel[1] = type_p[comp_idx - 1];
+ outLabel[2] = 0;
+ }
+ break;
+ case kGTOType_D:
+ outLabel[0] = 'D';
+ strncpy(outLabel + 1, type_d + comp_idx * 2, 2);
+ outLabel[3] = 0;
+ break;
+ case kGTOType_D5:
+ outLabel[0] = 'D';
+ strcpy(outLabel + 1, type_d5[comp_idx]);
+ break;
+ case kGTOType_F:
+ outLabel[0] = 'F';
+ strncpy(outLabel + 1, type_f + comp_idx * 3, 3);
+ outLabel[4] = 0;
+ break;
+ case kGTOType_F7:
+ outLabel[0] = 'F';
+ strcpy(outLabel + 1, type_f7[comp_idx]);
+ break;
+ case kGTOType_G:
+ outLabel[0] = 'G';
+ strcpy(outLabel + 1, type_g[comp_idx]);
+ break;
+ case kGTOType_G9:
+ outLabel[0] = 'G';
+ strcpy(outLabel + 1, type_g9[comp_idx]);
+ break;
+ default:
+ return -3; /* Unsupported orbital type (internal error) */
+ }
+ return 0;
+ }
+ }
+ return -4; /* comp_idx out of range? (internal error) */
+}
+
+/* Set MO coefficients for idx-th MO (1-based) */
int
MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Double *coeffs)
{
bset->nmos = bset->ncomps;
if (bset->nmos <= 0)
return -3; /* Bad or inconsistent number of MOs */
- bset->mo = (Double *)calloc(sizeof(Double), bset->nmos * bset->ncomps);
- bset->moenergies = (Double *)calloc(sizeof(Double), bset->nmos);
+ bset->mo = (Double *)calloc(sizeof(Double), (bset->nmos + 1) * bset->ncomps);
+ bset->moenergies = (Double *)calloc(sizeof(Double), bset->nmos + 1);
if (bset->mo == NULL || bset->moenergies == NULL) {
if (bset->mo != NULL)
free(bset->mo);
return -2; /* Low memory */
}
}
- if (idx < 0 || idx >= bset->nmos)
+ if (idx < 0)
+ idx = -idx + bset->ncomps;
+ if (idx < 0 || idx > bset->nmos)
return -4; /* Bad MO index */
+ if (idx == 0)
+ idx = bset->nmos; /* Arbitrary vector */
+ else
+ idx--;
if (energy != -1000000)
bset->moenergies[idx] = energy;
if (ncomps < bset->ncomps)
return 0;
}
-/* Allocate BasisSet record. rflag: UHF, 0; RHF, 1; ROHF, 2
- ne_alpha: number of alpha electrons, ne_beta: number of beta electrons
- The natoms and pos are copied from mol. */
+/* Get MO coefficients for idx-th MO (1-based) */
+/* Caution: *ncoeffs and *coeffs should be valid _before_ calling this function, i.e. */
+/* *ncoeffs = 0 && *coeffs = NULL or *coeffs is a valid memory pointer and *ncoeffs */
+/* properly designates the memory size as an array of Doubles. */
int
-MoleculeAllocateBasisSetRecord(Molecule *mol, Int rflag, Int ne_alpha, Int ne_beta)
+MoleculeGetMOCoefficients(Molecule *mol, Int idx, Double *energy, Int *ncoeffs, Double **coeffs)
+{
+ BasisSet *bset;
+ if (mol == NULL)
+ return -1; /* Molecule is empty */
+ bset = mol->bset;
+ if (bset == NULL || bset->ncomps <= 0)
+ return -2; /* No basis set info */
+ if (idx < 0)
+ idx = -idx + bset->ncomps;
+ if (idx < 0 || idx > bset->nmos)
+ return -3; /* MO index out of range */
+ if (idx == 0)
+ idx = bset->nmos; /* Arbitrary vector */
+ else
+ idx--;
+ if (energy != NULL)
+ *energy = bset->moenergies[idx];
+ if (ncoeffs != NULL && coeffs != NULL) {
+ if (*ncoeffs < bset->ncomps || *coeffs == NULL) {
+ if (*coeffs != NULL)
+ free(*coeffs); /* Caution: possible cause of SIGBUS if *coeff is not initialized properly */
+ *coeffs = (Double *)calloc(sizeof(Double), bset->ncomps);
+ *ncoeffs = bset->ncomps;
+ }
+ memmove(*coeffs, bset->mo + (idx * bset->ncomps), sizeof(Double) * bset->ncomps);
+ }
+ return 0;
+}
+
+/* Set Basic MO Info. rflag: 0, UHF; 1, RHF; 2, ROHF; -1, clear
+ ne_alpha: number of alpha electrons, ne_beta: number of beta electrons */
+int
+MoleculeSetMOInfo(Molecule *mol, Int rflag, Int ne_alpha, Int ne_beta)
{
BasisSet *bset;
- int i;
- Atom *ap;
if (mol == NULL || mol->natoms == 0)
return -1; /* Molecule is empty */
+ if (rflag < 0) {
+ if (mol->bset != NULL) {
+ BasisSetRelease(mol->bset);
+ mol->bset = NULL;
+ }
+ return 0;
+ }
bset = mol->bset;
if (bset == NULL) {
bset = mol->bset = (BasisSet *)calloc(sizeof(BasisSet), 1);
if (bset == NULL)
return -2; /* Low memory */
}
- if (bset->pos != NULL) {
- free(bset->pos);
- bset->pos = NULL;
- }
- bset->natoms = mol->natoms;
- bset->pos = (Vector *)calloc(sizeof(Vector), bset->natoms);
- if (bset->pos == NULL)
- return -2; /* Low memory */
- for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
- bset->pos[i].x = ap->r.x * kAngstrom2Bohr;
- bset->pos[i].y = ap->r.y * kAngstrom2Bohr;
- bset->pos[i].z = ap->r.z * kAngstrom2Bohr;
- }
+ bset->natoms_bs = mol->natoms;
bset->ne_alpha = ne_alpha;
bset->ne_beta = ne_beta;
bset->rflag = rflag;
dp[3] = d * 1.425410941;
dp += 5;
break;
+ /* TODO: Support F/F7 and G/G9 type orbitals */
}
}
}
Int *iary;
Double *dary;
Atom *ap;
- Vector *vp;
+/* Vector *vp; */
Double w;
*errbuf = NULL;
retval = 2;
goto cleanup;
}
+ bset->natoms_bs = natoms;
/* Allocate atom records (all are empty for now) */
AssignArray(&mp->atoms, &mp->natoms, gSizeOfAtomRecord, natoms - 1, NULL);
/* Also allocate atom position array for MO calculations */
- AssignArray(&bset->pos, &bset->natoms, sizeof(Vector), natoms - 1, NULL);
+ /* AssignArray(&bset->pos, &bset->natoms, sizeof(Vector), natoms - 1, NULL); */
/* Also allocate nuclear charge array */
bset->nuccharges = (Double *)calloc(sizeof(Double), natoms);
} else if (strcmp(buf, "Number of electrons") == 0) {
retval = 2;
goto cleanup;
}
- for (i = 0, ap = mp->atoms, vp = bset->pos; i < natoms; i++, ap = ATOM_NEXT(ap), vp++) {
- vp->x = dary[i * 3];
- vp->y = dary[i * 3 + 1];
- vp->z = dary[i * 3 + 2];
- ap->r.x = vp->x * kBohr2Angstrom;
- ap->r.y = vp->y * kBohr2Angstrom;
- ap->r.z = vp->z * kBohr2Angstrom;
+ for (i = 0, ap = mp->atoms; i < natoms; i++, ap = ATOM_NEXT(ap)) {
+ ap->r.x = dary[i * 3] * kBohr2Angstrom;
+ ap->r.y = dary[i * 3 + 1] * kBohr2Angstrom;
+ ap->r.z = dary[i * 3 + 2] * kBohr2Angstrom;
}
free(dary);
dary = NULL;
case -1: sp->sym = kGTOType_SP; sp->ncomp = 4; break;
case 2: sp->sym = kGTOType_D; sp->ncomp = 6; break;
case -2: sp->sym = kGTOType_D5; sp->ncomp = 5; break;
- /* TODO: Support F/F7 type orbitals */
- /* case 3: sp->sym = kGTOtype_F; sp->ncomp = 10; break;
- case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break; */
+ case 3: sp->sym = kGTOType_F; sp->ncomp = 10; break;
+ case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break;
+ case 4: sp->sym = kGTOType_G; sp->ncomp = 15; break;
+ case -4: sp->sym = kGTOType_G9; sp->ncomp = 9; break;
default:
s_append_asprintf(errbuf, "Line %d: unsupported shell type %d", lineNumber, iary[i]);
retval = 2;
vbuf[i].z = dval[3];
}
/* Skip until a blank line is found */
+ /* 2013.6.11. Line including "PM3" is also recognized as the end of atom */
while ((status = sReadLineWithInterrupt(buf, sizeof buf, fp, &lineNumber)) > 0) {
for (j = 0; buf[j] == ' '; j++);
- if (buf[j] == '\n')
+ if (buf[j] == '\n' || strncmp(buf + j, "PM3", 3) == 0)
break;
}
i++;
continue;
} else if (strstr(buf, "E(UHF)") != NULL || (strstr(buf, "E(RHF)") != NULL && (n1 = 1)) || (strstr(buf, "E(ROHF)") != NULL && (n1 = 2))) {
if (mol->bset == NULL) {
- i = MoleculeAllocateBasisSetRecord(mol, n1, 0, 0);
+ i = MoleculeSetMOInfo(mol, n1, 0, 0);
if (i != 0) {
s_append_asprintf(errbuf, "Line %d: cannot allocate basis set internal buffer", lineNumber);
retval = 8;
}
if (k < mol->bset->ncomps)
continue;
- j = MoleculeSetMOCoefficients(mol, i, -1000000, k, coeffs);
+ j = MoleculeSetMOCoefficients(mol, i + 1, -1000000, k, coeffs);
if (j != 0) {
s_append_asprintf(errbuf, "Line %d: cannot set coefficients for MO %d", lineNumber, i + 1);
free(coeffs);
if (newmol && mol->nbonds == 0) {
/* Guess bonds */
Int nbonds, *bonds;
- MoleculeGuessBonds(mol, 1.2, &nbonds, &bonds);
+ MoleculeGuessBonds(mol, 0.0, &nbonds, &bonds);
if (nbonds > 0) {
MolActionCreateAndPerform(mol, gMolActionAddBonds, nbonds * 2, bonds, NULL);
free(bonds);
MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
{
FILE *fp;
- Int i, j, k, n1, n2, n3, n_aniso, nframes;
+ Int i, j, k, n1, n2, n3, n_aniso, nframes, nanchors, n_uff;
Atom *ap;
+ char *p;
char bufs[6][8];
*errbuf = NULL;
fprintf(fp, "!:atoms\n");
fprintf(fp, "! idx seg_name res_seq res_name name type charge weight element atomic_number occupancy temp_factor int_charge\n");
- n1 = n2 = n3 = n_aniso = 0;
+ n1 = n2 = n3 = n_aniso = nanchors = n_uff = 0;
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
strncpy(bufs[0], ap->segName, 4);
bufs[0][4] = 0;
n3++;
if (ap->aniso != NULL)
n_aniso++;
+ if (ap->anchor != NULL)
+ nanchors++;
+ if (ap->uff_type[0] != 0)
+ n_uff++;
fprintf(fp, "%d %s %d %s %s %s %.5f %.5f %s %d %f %f %d\n", i, bufs[0], ap->resSeq, bufs[1], bufs[2], bufs[3], ap->charge, ap->weight, bufs[4], ap->atomicNumber, ap->occupancy, ap->tempFactor, ap->intCharge);
}
fprintf(fp, "\n");
+ if (n_uff > 0) {
+ fprintf(fp, "!:uff_type\n");
+ fprintf(fp, "! idx uff_type\n");
+ for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ fprintf(fp, "%d %.5s\n", i, ap->uff_type);
+ }
+ fprintf(fp, "\n");
+ }
+
if (n1 > 0) {
fprintf(fp, "!:atoms_symop\n");
fprintf(fp, "! idx symop symbase\n");
fprintf(fp, "\n");
}
+ if (nanchors > 0) {
+ fprintf(fp, "!:pi_anchor\n");
+ fprintf(fp, "! idx count; n1 weight1; n2 weight2; ...; nN weightN\n");
+ for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ Int *ip;
+ if (ap->anchor == NULL)
+ continue;
+ k = ap->anchor->connect.count;
+ ip = AtomConnectData(&ap->anchor->connect);
+ fprintf(fp, "%d %d\n", i, k);
+ for (j = 0; j < k; j++) {
+ fprintf(fp, "%d %f\n", ip[j], ap->anchor->coeffs[j]);
+ }
+ }
+ fprintf(fp, "\n");
+ }
+
n1 = nframes;
if (n1 > 0)
n2 = mp->cframe;
}
fprintf(fp, "\n");
}
+
+ if (mp->nbondOrders > 0) {
+ fprintf(fp, "!:bond_orders\n");
+ fprintf(fp, "! order1 order2 order3 order4\n");
+ for (i = 0; i < mp->nbondOrders; i++) {
+ fprintf(fp, "%.6f%c", mp->bondOrders[i], (i % 4 == 3 || i == mp->nbondOrders - 1 ? '\n' : ' '));
+ }
+ fprintf(fp, "\n");
+ }
if (mp->nangles > 0) {
fprintf(fp, "!:angles\n");
fprintf(fp, "\n");
}
-#if PIATOM
- if (mp->npiatoms > 0) {
- PiAtom *pp;
- fprintf(fp, "!:pi_atoms\n");
- fprintf(fp, "! name type n; a1 coeff1; a2 coeff2; ...; a_n coeff_n\n");
- for (i = 0, pp = mp->piatoms; i < mp->npiatoms; i++, pp++) {
- strncpy(bufs[0], pp->aname, 4);
- bufs[0][4] = 0;
- AtomTypeDecodeToString(pp->type, bufs[1]);
- bufs[1][6] = 0;
- for (j = 0; j < 2; j++) {
- if (bufs[j][0] == 0) {
- bufs[j][0] = '_';
- bufs[j][1] = 0;
- }
- }
- fprintf(fp, "%s %s %d\n", bufs[0], bufs[1], pp->connect.count);
- ip = AtomConnectData(&pp->connect);
- for (j = 0; j < pp->connect.count; j++) {
- fprintf(fp, "%d %g\n", ip[j], (j < pp->ncoeffs ? pp->coeffs[j] : 0.0));
- }
- }
- fprintf(fp, "\n");
- }
-
- if (mp->npibonds > 0) {
- fprintf(fp, "!:pi_atom_constructs\n");
- fprintf(fp, "! a1 b1 c1 d1 a2 b2 c2 d2\n");
- for (i = 0; i < mp->npibonds * 4; i++) {
- j = mp->pibonds[i];
- if (j >= ATOMS_MAX_NUMBER)
- j = -2 - (j - ATOMS_MAX_NUMBER);
- else if (j < 0)
- j = -1;
- fprintf(fp, "%d%c", j, (i % 8 == 7 || i == mp->npibonds * 4 - 1 ? '\n' : ' '));
- }
- fprintf(fp, "\n");
- }
-#endif /* PIATOM */
-
if (n_aniso > 0) {
fprintf(fp, "!:anisotropic_thermal_parameters\n");
fprintf(fp, "! b11 b22 b33 b12 b13 b23 [sigma; sb11 sb22 sb33 sb12 sb13 sb23]\n");
fprintf(fp, "cutoff %g\n", arena->cutoff);
fprintf(fp, "electro_cutoff %g\n", arena->electro_cutoff);
fprintf(fp, "pairlist_distance %g\n", arena->pairlist_distance);
+ fprintf(fp, "switch_distance %g\n", arena->switch_distance);
fprintf(fp, "temperature %g\n", arena->temperature);
fprintf(fp, "andersen_freq %d\n", arena->andersen_thermo_freq);
fprintf(fp, "andersen_coupling %g\n", arena->andersen_thermo_coupling);
}
if (mp->mview != NULL) {
- float f[4];
+ double f[4];
if (mp->mview->track != NULL) {
fprintf(fp, "!:trackball\n");
fprintf(fp, "! scale; trx try trz; theta_deg x y z\n");
mp->mview->showPeriodicImage[0], mp->mview->showPeriodicImage[1],
mp->mview->showPeriodicImage[2], mp->mview->showPeriodicImage[3],
mp->mview->showPeriodicImage[4], mp->mview->showPeriodicImage[5]);
+ if (mp->mview->atomRadius != 0.2)
+ fprintf(fp, "atom_radius %f\n", mp->mview->atomRadius);
+ if (mp->mview->bondRadius != 0.1)
+ fprintf(fp, "bond_radius %f\n", mp->mview->bondRadius);
+ if (mp->mview->atomResolution != 12)
+ fprintf(fp, "atom_resolution %d\n", mp->mview->atomResolution);
+ if (mp->mview->bondResolution != 8)
+ fprintf(fp, "bond_resolution %d\n", mp->mview->bondResolution);
+ fprintf(fp, "\n");
+ }
+
+ if (mp->nmolprops > 0) {
+ MolProp *prp;
+ for (i = 0, prp = mp->molprops; i < mp->nmolprops; i++, prp++) {
+ /* Encode the property name if necessary */
+ char enc[1024];
+ n1 = n2 = 0;
+ for (p = prp->propname; *p != 0 && n1 < 900; p++) {
+ if (*p > ' ' && *p != '%' && *p < 0x7f) {
+ enc[n1++] = *p;
+ n2 = n1;
+ } else {
+ sprintf(enc + n1, "%%%02x", *p);
+ n1 += 3;
+ }
+ }
+ if (*p == 0)
+ enc[n1] = 0;
+ else {
+ enc[n2] = 0; /* Truncate after last ASCII character */
+ n1 = n2;
+ }
+ if (n1 == 0) {
+ sprintf(enc, "prop_%d", i + 1);
+ n1 = strlen(enc);
+ }
+ fprintf(fp, "!:property ; %s\n", enc);
+ for (j = 0; j < nframes; j++) {
+ fprintf(fp, "%.18g\n", prp->propvals[j]);
+ }
+ fprintf(fp, "\n");
+ }
+ }
+
+ if (mp->bset != NULL) {
+ /* Gaussian primitive info */
+ ShellInfo *sp;
+ PrimInfo *pp;
+ fprintf(fp, "!:gaussian_primitives\n");
+ fprintf(fp, "! sym nprims a_idx; A C Csp\n");
+ for (i = 0, sp = mp->bset->shells; i < mp->bset->nshells; i++, sp++) {
+ switch (sp->sym) {
+ case kGTOType_S: p = "S"; break;
+ case kGTOType_P: p = "P"; break;
+ case kGTOType_SP: p = "SP"; break;
+ case kGTOType_D: p = "D"; break;
+ case kGTOType_D5: p = "D5"; break;
+ case kGTOType_F: p = "F"; break;
+ case kGTOType_F7: p = "F7"; break;
+ case kGTOType_G: p = "G"; break;
+ case kGTOType_G9: p = "G9"; break;
+ default: snprintf(bufs[0], 8, "X%d", sp->sym); p = bufs[0]; break;
+ }
+ fprintf(fp, "%s %d %d\n", p, sp->nprim, sp->a_idx);
+ pp = mp->bset->priminfos + sp->p_idx;
+ for (j = 0; j < sp->nprim; j++, pp++) {
+ fprintf(fp, "%.18g %.18g %.18g\n", pp->A, pp->C, pp->Csp);
+ }
+ }
+ fprintf(fp, "\n");
+
+ /* MO info */
+ fprintf(fp, "!:mo_info\n");
+ fprintf(fp, "! uhf|rhf|rohf ne_alpha ne_beta\n");
+ switch (mp->bset->rflag) {
+ case 0: p = "UHF"; break;
+ case 1: p = "RHF"; break;
+ case 2: p = "ROHF"; break;
+ default: p = "(unknown)"; break;
+ }
+ fprintf(fp, "%s %d %d\n", p, mp->bset->ne_alpha, mp->bset->ne_beta);
+ fprintf(fp, "\n");
+
+ /* MO coefficients */
+ fprintf(fp, "!:mo_coefficients\n");
+ for (i = 0; i < mp->bset->nmos; i++) {
+ fprintf(fp, "MO %d %.18g\n", i + 1, mp->bset->moenergies[i]);
+ for (j = 0; j < mp->bset->ncomps; j++) {
+ fprintf(fp, "%.18g%c", mp->bset->mo[i * mp->bset->ncomps + j], (j % 6 == 5 || j == mp->bset->ncomps - 1 ? '\n' : ' '));
+ }
+ }
fprintf(fp, "\n");
}
+ if (mp->mview != NULL && mp->mview->ngraphics > 0) {
+ MainViewGraphic *gp;
+ fprintf(fp, "!:graphics\n");
+ for (i = 0; i < mp->mview->ngraphics; i++) {
+ gp = mp->mview->graphics + i;
+ switch (gp->kind) {
+ case kMainViewGraphicLine: fprintf(fp, "line\n"); break;
+ case kMainViewGraphicPoly: fprintf(fp, "poly\n"); break;
+ case kMainViewGraphicCylinder: fprintf(fp, "cylinder\n"); break;
+ case kMainViewGraphicCone: fprintf(fp, "cone\n"); break;
+ case kMainViewGraphicEllipsoid: fprintf(fp, "ellipsoid\n"); break;
+ default: fprintf(fp, "unknown\n"); break;
+ }
+ fprintf(fp, "%d %d\n", gp->closed, gp->visible);
+ fprintf(fp, "%.4f %.4f %.4f %.4f\n", gp->rgba[0], gp->rgba[1], gp->rgba[2], gp->rgba[3]);
+ fprintf(fp, "%d\n", gp->npoints);
+ for (j = 0; j < gp->npoints; j++)
+ fprintf(fp, "%.6f %.6f %.6f\n", gp->points[j * 3], gp->points[j * 3 + 1], gp->points[j * 3 + 2]);
+ fprintf(fp, "%d\n", gp->nnormals);
+ for (j = 0; j < gp->nnormals; j++)
+ fprintf(fp, "%.6f %.6f %.6f\n", gp->normals[j * 3], gp->normals[j * 3 + 1], gp->normals[j * 3 + 2]);
+ }
+ fprintf(fp, "\n");
+ }
+
+ /* Plug-in in the Ruby world */
+ {
+ char *outMessage;
+ if (MolActionCreateAndPerform(mp, SCRIPT_ACTION(";s"),
+ "proc { savembsf_plugin rescue \"Plug-in error: #{$!.to_s}\" }", &outMessage) == 0) {
+ if (outMessage[0] != 0) {
+ if (strncmp(outMessage, "Plug-in", 7) == 0) {
+ s_append_asprintf(errbuf, "%s", outMessage);
+ } else {
+ fprintf(fp, "%s\n", outMessage);
+ }
+ }
+ free(outMessage);
+ }
+ }
+
fclose(fp);
return 0;
}
fprintf(fp, " %.8g %.8g %.8g ! %d,%.4s\n", r.x, r.y, r.z, i + 1, ap->aname);
}
fprintf(fp, "\n");
-#if 0
- if (mp->nframes > 0) {
- int fn; /* Frame number */
- for (fn = 0; fn < ap->nframes; fn++) {
- fprintf(fp, "%8d !COORD: coordinates for frame %d\n", mp->natoms, fn);
- for (i = 0; i < mp->natoms; i++) {
- Vector r;
- ap = ATOM_AT_INDEX(mp->atoms, i);
- if (ap->frames == NULL || fn >= ap->nframes)
- r = ap->r;
- else
- r = ap->frames[fn];
- fprintf(fp, " %.8g %.8g %.8g ! %d,%.4s\n", r.x, r.y, r.z, i + 1, ap->name);
- }
- fprintf(fp, "\n");
- }
- }
-#endif
}
fclose(fp);
free(exbonds);
}
}
-
-#if 0
-{
- /* Explicit bond table, sorted by bond type */
- for (i = j = 0; i < mp->nbonds; i++) {
- n1 = mp->bonds[i * 2];
- n2 = mp->bonds[i * 2 + 1];
- ap1 = ATOM_AT_INDEX(mp->atoms, n1);
- ap2 = ATOM_AT_INDEX(mp->atoms, n2);
- if ((ap1->exflags & kAtomHiddenFlag) || (ap2->exflags & kAtomHiddenFlag))
- continue;
- if (ap1->atomicNumber > 18 || ap2->atomicNumber > 18) {
- type = 3;
- } else if (ap1->atomicNumber > 1 && ap1->atomicNumber > 1) {
- type = 2;
- } else {
- type = 1;
- }
- ip[j * 3] = type;
- ip[j * 3 + 1] = sMakeAdc(n1, ap1->symbase, ap1->symop);
- ip[j * 3 + 2] = sMakeAdc(n2, ap2->symbase, ap2->symop);
- j++;
- }
- mergesort(ip, j, sizeof(int) * 3, sCompareBondType);
-
- /* Output instruction cards */
- strcpy(buf, " 1 811");
- for (i = n1 = 0; i < j; i++) {
- n2 = (n1 % 3) * 18 + 9;
- snprintf(buf + n2, 80 - n2, "%9d%9d\n", ip[i * 3 + 1], ip[i * 3 + 2]);
- if (i == j - 1 || n1 >= 29 || ip[i * 3] != ip[i * 3 + 3]) {
- /* End of this instruction */
- buf[2] = '2';
- fputs(buf, fp);
- switch (ip[i * 3]) {
- case 3: rad = 0.06; nshades = 5; break;
- case 2: rad = 0.06; nshades = 1; break;
- default: rad = 0.04; nshades = 1; break;
- }
- fprintf(fp, " %3d %6.3f\n", nshades, rad);
- strcpy(buf, " 1 811");
- n1 = 0;
- continue;
- } else if (n1 % 3 == 2) {
- fputs(buf, fp);
- strcpy(buf, " 1 ");
- }
- n1++;
- }
- free(ip);
-}
-#endif
int
MoleculeWriteToTepFile(Molecule *mp, const char *fname, char **errbuf)
IntGroupRelease(ig3);
}
+ {
+ /* Update the path information of the molecule before MD setup */
+ char *buf = (char *)malloc(4096);
+ MoleculeCallback_pathName(mol, buf, sizeof buf);
+ MoleculeSetPath(mol, buf);
+ free(buf);
+ }
+
/* Prepare parameters and internal information */
msg = md_prepare(arena, check_only);
{
Molecule *mp;
Parameter *par;
+ Atom *ap;
/* int result; */
mp = MoleculeNew();
NewArray(&mp->atoms, &mp->natoms, gSizeOfAtomRecord, n);
memmove(mp->atoms, ptr, len);
} else if (strcmp(data, "ANISO") == 0) {
- Atom *ap;
n = len / (sizeof(Int) + sizeof(Aniso));
for (i = 0; i < n; i++) {
j = *((const Int *)ptr);
ptr += sizeof(Int) + sizeof(Aniso);
}
} else if (strcmp(data, "FRAME") == 0) {
- Atom *ap;
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
if (ap->nframes == 0)
continue;
- ap->frames = (Vector *)malloc(sizeof(Vector) * ap->nframes);
+ n = ap->nframes;
+ ap->frames = NULL;
+ ap->nframes = 0;
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), n);
if (ap->frames == NULL)
goto out_of_memory;
memmove(ap->frames, ptr, sizeof(Vector) * ap->nframes);
ptr += sizeof(Vector) * ap->nframes;
}
} else if (strcmp(data, "EXTCON") == 0) {
- Atom *ap;
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
if (ap->connect.count <= ATOM_CONNECT_LIMIT)
continue;
n = len / sizeof(Transform);
NewArray(&mp->syms, &mp->nsyms, sizeof(Transform), n);
memmove(mp->syms, ptr, len);
-#if PIATOM
- } else if (strcmp(data, "PIATOM") == 0) {
+ } else if (strcmp(data, "ANCHOR") == 0) {
const char *ptr2 = ptr + len;
- mp->piatoms = NULL;
- mp->npiatoms = 0;
while (ptr < ptr2) {
- PiAtom pa;
- memset(&pa, 0, sizeof(pa));
- n = offsetof(PiAtom, connect);
- memmove(&pa, ptr, n);
- ptr += n;
- n = *((Int *)ptr);
- if (n > 0) {
- AtomConnectResize(&pa.connect, n);
- memmove(AtomConnectData(&pa.connect), ptr + sizeof(Int), n * sizeof(Int));
- }
- ptr += sizeof(Int) * (n + 1);
- n = *((Int *)ptr);
- if (n > 0) {
- NewArray(&pa.coeffs, &pa.ncoeffs, sizeof(Double), n);
- memmove(pa.coeffs, ptr + sizeof(Int), n * sizeof(Double));
- }
- ptr += sizeof(Int) + sizeof(Double) * n;
- AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, &pa);
- }
- } else if (strcmp(data, "PIBOND") == 0) {
- n = len / (sizeof(Int) * 4);
- NewArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, n);
- memmove(mp->pibonds, ptr, len);
-#endif /* PIATOM */
+ PiAnchor an;
+ memset(&an, 0, sizeof(an));
+ i = *((Int *)ptr);
+ if (i >= 0 && i < mp->natoms) {
+ n = *((Int *)(ptr + sizeof(Int)));
+ AtomConnectResize(&(an.connect), n);
+ memmove(AtomConnectData(&(an.connect)), ptr + sizeof(Int) * 2, sizeof(Int) * n);
+ NewArray(&an.coeffs, &an.ncoeffs, sizeof(Double), n);
+ memmove(an.coeffs, ptr + sizeof(Int) * (2 + n), sizeof(Double) * n);
+ ap = ATOM_AT_INDEX(mp->atoms, i);
+ ap->anchor = (PiAnchor *)malloc(sizeof(PiAnchor));
+ memmove(ap->anchor, &an, sizeof(PiAnchor));
+ }
+ ptr += sizeof(Int) * (2 + n) + sizeof(Double) * n;
+ }
} else if (strcmp(data, "TIME") == 0) {
if (timep != NULL)
*timep = *((Int *)ptr);
MoleculeSerialize(Molecule *mp, Int *outLength, Int *timep)
{
char *ptr, *p;
- int len, len_all, i, naniso, nframes, nconnects;
+ int len, len_all, i, naniso, nframes, nconnects, nanchors;
+ Atom *ap;
/* Array of atoms */
len = 8 + sizeof(Int) + gSizeOfAtomRecord * mp->natoms;
*((Int *)(ptr + 8)) = gSizeOfAtomRecord * mp->natoms;
p = ptr + 8 + sizeof(Int);
memmove(p, mp->atoms, gSizeOfAtomRecord * mp->natoms);
- naniso = nframes = nconnects = 0;
+ naniso = nframes = nconnects = nanchors = 0;
for (i = 0; i < mp->natoms; i++) {
- Atom *ap = ATOM_AT_INDEX(p, i);
+ ap = ATOM_AT_INDEX(p, i);
if (ap->aniso != NULL) {
naniso++;
ap->aniso = NULL;
nconnects += ap->connect.count;
ap->connect.u.ptr = NULL;
}
+ if (ap->anchor != NULL) {
+ nanchors++;
+ ap->anchor = NULL;
+ }
}
len_all = len;
*((Int *)(p + 8)) = (sizeof(Int) + sizeof(Aniso)) * naniso;
p += 8 + sizeof(Int);
for (i = 0; i < mp->natoms; i++) {
- Atom *ap = ATOM_AT_INDEX(mp->atoms, i);
+ ap = ATOM_AT_INDEX(mp->atoms, i);
if (ap->aniso != NULL) {
*((Int *)p) = i;
*((Aniso *)(p + sizeof(Int))) = *(ap->aniso);
*((Int *)(p + 8)) = sizeof(Vector) * nframes;
p += 8 + sizeof(Int);
for (i = 0; i < mp->natoms; i++) {
- Atom *ap = ATOM_AT_INDEX(mp->atoms, i);
+ ap = ATOM_AT_INDEX(mp->atoms, i);
if (ap->frames != NULL) {
memmove(p, ap->frames, sizeof(Vector) * ap->nframes);
p += sizeof(Vector) * ap->nframes;
*((Int *)(p + 8)) = sizeof(Int) * nconnects;
p += 8 + sizeof(Int);
for (i = 0; i < mp->natoms; i++) {
- Atom *ap = ATOM_AT_INDEX(mp->atoms, i);
+ ap = ATOM_AT_INDEX(mp->atoms, i);
if (ap->connect.count > ATOM_CONNECT_LIMIT) {
memmove(p, ap->connect.u.ptr, sizeof(Int) * ap->connect.count);
p += sizeof(Int) * ap->connect.count;
len_all += len;
}
-#if PIATOM
- /* Pi-atoms */
- if (mp->npiatoms > 0) {
+ /* Pi-anchors */
+ if (nanchors > 0) {
/* Estimate the necessary storage first */
+ /* One entry consists of { atom_index (Int), number_of_connects (Int), connects (Int's), weights (Double's) } */
len = 8 + sizeof(Int);
- for (i = 0; i < mp->npiatoms; i++) {
- len += offsetof(PiAtom, connect); /* Members before 'connect' is stored as they are */
- len += sizeof(Int) * (1 + mp->piatoms[i].connect.count); /* Array of Int's */
- len += sizeof(Int) + sizeof(Double) * mp->piatoms[i].ncoeffs; /* Array of Double's */
+ for (i = 0; i < mp->natoms; i++) {
+ ap = ATOM_AT_INDEX(mp->atoms, i);
+ if (ap->anchor != NULL)
+ len += sizeof(Int) * 2 + (sizeof(Int) + sizeof(Double)) * ap->anchor->connect.count;
}
ptr = (char *)realloc(ptr, len_all + len);
if (ptr == NULL)
goto out_of_memory;
p = ptr + len_all;
- memmove(p, "PIATOM\0\0", 8);
+ memmove(p, "ANCHOR\0\0", 8);
*((Int *)(p + 8)) = len - (8 + sizeof(Int));
p += 8 + sizeof(Int);
- for (i = 0; i < mp->npiatoms; i++) {
- int len0;
- PiAtom *pp = &(mp->piatoms[i]);
- len0 = offsetof(PiAtom, connect);
- memmove(p, pp, len0);
- p += len0;
- len0 = pp->connect.count * sizeof(Int);
- *((Int *)p) = pp->connect.count;
- if (len0 > 0)
- memmove(p + sizeof(Int), AtomConnectData(&(pp->connect)), len0);
- p += sizeof(Int) + len0;
- len0 = pp->ncoeffs * sizeof(Double);
- *((Int *)p) = pp->ncoeffs;
- if (len0 > 0)
- memmove(p + sizeof(Int), pp->coeffs, len0);
- p += sizeof(Int) + len0;
+ for (i = 0; i < mp->natoms; i++) {
+ Int count, *ip;
+ ap = ATOM_AT_INDEX(mp->atoms, i);
+ if (ap->anchor != NULL) {
+ count = ap->anchor->connect.count;
+ *((Int *)p) = i;
+ *((Int *)(p + sizeof(Int))) = count;
+ p += sizeof(Int) * 2;
+ ip = AtomConnectData(&(ap->anchor->connect));
+ memmove(p, ip, sizeof(Int) * count);
+ p += sizeof(Int) * count;
+ memmove(p, ap->anchor->coeffs, sizeof(Double) * count);
+ p += sizeof(Double) * count;
+ }
}
len_all += len;
}
- /* Pi-atom constructs */
- if (mp->npibonds > 0) {
- len = 8 + sizeof(Int) + sizeof(Int) * 4 * mp->npibonds;
- ptr = (char *)realloc(ptr, len_all + len);
- if (ptr == NULL)
- goto out_of_memory;
- p = ptr + len_all;
- memmove(p, "PIBOND\0\0", 8);
- *((Int *)(p + 8)) = sizeof(Int) * 4 * mp->npibonds;
- p += 8 + sizeof(Int);
- memmove(p, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
- len_all += len;
- }
-#endif /* PIATOM */
-
/* Parameters */
if (mp->par != NULL) {
int type;
return sMoleculeSearchAcrossAtomGroup(mp->nimpropers, mp->impropers, 4, atomgroup, "impropers");
}
-/* Subroutine for MoleculeGuessBonds. It can be also used independently, but make sure that *outNbonds/*outBonds
+/* 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);
- }
- }
- */
+ if (limit == 0.0)
+ limit = 1.2;
+ 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;
return retval;
}
+int
+MoleculeAreAtomsConnected(Molecule *mol, int idx1, int idx2)
+{
+ Atom *ap1 = ATOM_AT_INDEX(mol->atoms, idx1);
+ if (AtomConnectHasEntry(&ap1->connect, idx2))
+ return 1;
+ else if (ap1->anchor != NULL && AtomConnectHasEntry(&(ap1->anchor->connect), idx2))
+ return 2;
+ else return 0;
+}
+
#pragma mark ====== Atom names ======
/* Look for the n1-th atom in resno-th residue (n1 is 0-based) */
return -1; /* Not found */
}
-int
-MoleculeAreAtomsConnected(Molecule *mp, int n1, int n2)
-{
- Atom *ap;
- Int i, *cp;
- if (mp == NULL || n1 < 0 || n1 >= mp->natoms || n2 < 0 || n2 >= mp->natoms)
- return 0;
- ap = ATOM_AT_INDEX(mp->atoms, n1);
- cp = AtomConnectData(&ap->connect);
- for (i = 0; i < ap->connect.count; i++)
- if (cp[i] == n2)
- return 1;
- return 0;
-}
-
-
void
MoleculeGetAtomName(Molecule *mp, int index, char *buf, int bufsize)
{
int i, j, n[3];
if (mp == NULL || mp->cell == NULL)
return -1;
- if (mp->nsyms == 0)
- return -2;
if (is_cartesian) {
TransformMul(t, tf, mp->cell->tr);
TransformMul(t, mp->cell->rtr, t);
} else {
memmove(t, tf, sizeof(Transform));
}
- for (i = 0; i < mp->nsyms; i++) {
- Transform *tp = mp->syms + i;
+ for (i = 0; i < mp->nsyms || i == 0; i++) {
+ Transform *tp = &(SYMMETRY_AT_INDEX(mp->syms, i));
for (j = 0; j < 9; j++) {
if (fabs((*tp)[j] - t[j]) > 1e-4)
break;
If indices is non-NULL, it should be an array of Int with at least
IntGroupGetCount(group) entries, and on return it contains the
indices of the expanded atoms (may be existing atoms if the expanded
- atoms are already present) */
+ atoms are already present)
+ If allowOverlap is non-zero, then the new atom is created even when the
+ coordinates coincide with the some other atom (special position) of the
+ same element; otherwise, such atom will not be created and the existing
+ atom is returned in indices[]. */
int
-MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indices)
+MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indices, Int allowOverlap)
{
- int i, n, n0, n1, count, *table;
+ int i, n, n0, n1, n2, base, count, *table;
Atom *ap;
IntGroupIterator iter;
- Transform tr;
-
+ Transform tr, t1;
+ Symop symop1;
+ Atom *ap2;
+ Vector nr, dr;
+
if (mp == NULL || mp->natoms == 0 || group == NULL || (count = IntGroupGetCount(group)) == 0)
return -1;
- if (symop.sym >= mp->nsyms)
+ if (symop.sym != 0 && symop.sym >= mp->nsyms)
return -2;
/* Create atoms, with avoiding duplicates */
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(mp);
-#endif
n0 = n1 = mp->natoms;
table = (int *)malloc(sizeof(int) * n0);
if (table == NULL)
MoleculeGetTransformForSymop(mp, symop, &tr, 0);
__MoleculeLock(mp);
for (i = 0; i < count; i++) {
- int n2, base;
- Symop symop1;
- Atom *ap2;
- Vector nr, dr;
n = IntGroupIteratorNext(&iter);
ap = ATOM_AT_INDEX(mp->atoms, n);
if (SYMOP_ALIVE(ap->symop)) {
/* Calculate the cumulative symop */
- Transform t1;
+ Transform tr2;
MoleculeGetTransformForSymop(mp, ap->symop, &t1, 0);
- TransformMul(t1, tr, t1);
- if (MoleculeGetSymopForTransform(mp, t1, &symop1, 0) != 0) {
+ TransformMul(tr2, tr, t1);
+ if (MoleculeGetSymopForTransform(mp, tr2, &symop1, 0) != 0) {
if (indices != NULL)
indices[i] = -1;
continue; /* Skip this atom */
symop1 = symop;
base = n;
}
+
+ /* Calculate the expande position */
+ MoleculeTransformBySymop(mp, &(ap->r), &nr, symop);
+
/* Is this expansion already present? */
for (n2 = 0, ap2 = mp->atoms; n2 < n0; n2++, ap2 = ATOM_NEXT(ap2)) {
+ /* Symmetry operation and the base atom are the same */
if (ap2->symbase == base && SYMOP_EQUAL(symop1, ap2->symop))
break;
+ /* Atomic number and the position are the same */
+ if (ap2->atomicNumber == ap->atomicNumber && allowOverlap == 0) {
+ VecSub(dr, ap2->r, nr);
+ if (VecLength2(dr) < 1e-6)
+ break;
+ }
}
if (n2 < n0) {
/* If yes, then skip it */
if (indices != NULL)
indices[i] = n2;
continue;
- }
- /* Is the expanded position coincides with itself? */
- MoleculeTransformBySymop(mp, &(ap->r), &nr, symop);
- VecSub(dr, ap->r, nr);
- if (VecLength2(dr) < 1e-6) {
- /* If yes, then this atom is included but no new atom is created */
- table[n] = n;
- if (indices != NULL)
- indices[i] = n;
} else {
/* Create a new atom */
Atom newAtom;
AtomClean(&newAtom);
ap2 = ATOM_AT_INDEX(mp->atoms, mp->natoms - 1);
ap2->r = nr;
- ap2->symbase = n;
- ap2->symop = symop;
- ap2->symop.alive = (symop.dx != 0 || symop.dy != 0 || symop.dz != 0 || symop.sym != 0);
+ ap2->symbase = base;
+ ap2->symop = symop1;
+ ap2->symop.alive = (symop1.dx != 0 || symop1.dy != 0 || symop1.dz != 0 || symop1.sym != 0);
table[n] = n1; /* The index of the new atom */
MoleculeSetAnisoBySymop(mp, n1); /* Recalculate anisotropic parameters according to symop */
if (indices != NULL)
IntGroupIteratorRelease(&iter);
/* Create bonds */
- for (i = 0; i < n0; i++) {
- int b[2];
- Int *cp;
- b[0] = table[i];
- if (b[0] < 0 || b[0] == i)
- continue;
+ for (i = n0; i < n1; i++) {
+ Int b[2], j;
ap = ATOM_AT_INDEX(mp->atoms, i);
- cp = AtomConnectData(&ap->connect);
- for (n = 0; n < ap->connect.count; n++) {
- b[1] = table[cp[n]];
- if (b[1] < 0)
- continue;
- if (b[1] > n0 && b[0] > b[1])
- continue;
- MoleculeAddBonds(mp, 1, b, NULL, 1);
+ if (SYMOP_ALIVE(ap->symop) && MoleculeGetTransformForSymop(mp, ap->symop, &tr, 1) == 0) {
+ /* For each connected atom, look for the transformed atom */
+ Int *cp;
+ ap2 = ATOM_AT_INDEX(mp->atoms, ap->symbase);
+ cp = AtomConnectData(&ap2->connect);
+ n2 = ap2->connect.count;
+ for (n = 0; n < n2; n++) {
+ Atom *apn = ATOM_AT_INDEX(mp->atoms, cp[n]);
+ nr = apn->r;
+ TransformVec(&nr, tr, &nr);
+ /* Look for the bonded atom transformed by ap->symop */
+ for (j = 0, ap2 = mp->atoms; j < mp->natoms; j++, ap2 = ATOM_NEXT(ap2)) {
+ if (ap2->symbase == cp[n] && SYMOP_EQUAL(ap->symop, ap2->symop))
+ break;
+ VecSub(dr, nr, ap2->r);
+ if (ap2->atomicNumber == apn->atomicNumber && VecLength2(dr) < 1e-6)
+ break;
+ }
+ if (j < mp->natoms) {
+ /* Bond i-j is created */
+ b[0] = i;
+ b[1] = j;
+ if (MoleculeLookupBond(mp, b[0], b[1]) < 0)
+ MoleculeAddBonds(mp, 1, b, NULL, 1);
+ }
+ }
}
}
mp->needsMDRebuild = 1;
}
/* Recalculate the coordinates of symmetry expanded atoms.
+ (Also recalculate the positions of pi-anchor atoms)
Returns the number of affected atoms.
If group is non-NULL, only the expanded atoms whose base atoms are in the
given group are considered.
{
int i, count;
Atom *ap, *bp;
- if (mp == NULL || mp->natoms == 0 || mp->nsyms == 0)
- return 0;
- if (groupout != NULL && vpout != NULL) {
- *groupout = IntGroupNew();
- if (*groupout == NULL)
- return -1;
- *vpout = (Vector *)malloc(sizeof(Vector) * mp->natoms);
- if (*vpout == NULL) {
- IntGroupRelease(*groupout);
- return -1;
- }
- } else groupout = NULL; /* To simplify test for validity of groupout/vpout */
+ Vector nr, dr;
+ IntGroup *ig = NULL;
+ Vector *vp = NULL;
+ if (mp == NULL || mp->natoms == 0)
+ return 0;
+
__MoleculeLock(mp);
count = 0;
+ if (mp->nsyms != 0) {
+ for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ if (!SYMOP_ALIVE(ap->symop))
+ continue;
+ if (group != NULL && IntGroupLookup(group, ap->symbase, NULL) == 0)
+ continue;
+ bp = ATOM_AT_INDEX(mp->atoms, ap->symbase);
+ MoleculeTransformBySymop(mp, &(bp->r), &nr, ap->symop);
+ VecSub(dr, nr, ap->r);
+ if (VecLength2(dr) < 1e-20)
+ continue;
+ if (groupout != NULL) {
+ if (ig == NULL) {
+ ig = IntGroupNew();
+ vp = (Vector *)calloc(sizeof(Vector), mp->natoms);
+ }
+ vp[count] = ap->r;
+ IntGroupAdd(ig, i, 1);
+ }
+ ap->r = nr;
+ count++;
+ }
+ }
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
- Vector nr, dr;
- if (!SYMOP_ALIVE(ap->symop))
+ Int *ip, j, n;
+ if (ap->anchor == NULL)
continue;
- if (group != NULL && IntGroupLookup(group, ap->symbase, NULL) == 0)
- continue;
- bp = ATOM_AT_INDEX(mp->atoms, ap->symbase);
- MoleculeTransformBySymop(mp, &(bp->r), &nr, ap->symop);
+ if (group != NULL) {
+ if (IntGroupLookup(group, i, NULL) == 0) {
+ n = ap->anchor->connect.count;
+ ip = AtomConnectData(&(ap->anchor->connect));
+ for (j = 0; j < n; j++) {
+ if (IntGroupLookup(group, ip[j], NULL) != 0)
+ break;
+ }
+ if (j == n)
+ continue; /* This pi-anchor should not be modified */
+ }
+ }
+ nr = ap->r;
+ MoleculeCalculatePiAnchorPosition(mp, i);
VecSub(dr, nr, ap->r);
- if (VecLength2(dr) < 1e-20)
+ if (VecLength2(dr) < 1e-20) {
+ ap->r = nr; /* No change */
continue;
+ }
if (groupout != NULL) {
- (*vpout)[count] = ap->r;
- IntGroupAdd(*groupout, i, 1);
+ if (ig == NULL) {
+ ig = IntGroupNew();
+ vp = (Vector *)calloc(sizeof(Vector), mp->natoms);
+ }
+ vp[count] = nr;
+ IntGroupAdd(ig, i, 1);
}
- ap->r = nr;
count++;
}
mp->needsMDCopyCoordinates = 1;
__MoleculeUnlock(mp);
- if (groupout != NULL) {
- if (count == 0) {
- free(*vpout);
- *vpout = NULL;
- IntGroupRelease(*groupout);
- *groupout = NULL;
+
+ if (count > 0) {
+ if (groupout != NULL && vpout != NULL) {
+ *groupout = ig;
+ *vpout = (Vector *)realloc(vp, sizeof(Vector) * count);
} else {
- *vpout = (Vector *)realloc(*vpout, sizeof(Vector) * count);
+ IntGroupRelease(ig);
+ free(vp);
+ }
+ } else {
+ if (groupout != NULL && vpout != NULL) {
+ *groupout = NULL;
+ *vpout = NULL;
}
}
return count;
sRemoveElementsFromArrayAtPositions(void *objs, int nobjs, void *clip, size_t size, IntGroup *where)
{
int n1, n2, n3, start, end, i;
- if (objs == NULL || where == NULL)
+ if (where == NULL || IntGroupGetCount(where) == 0)
+ return 0; /* No operation */
+ if (objs == NULL || nobjs == 0)
return 1; /* Bad argument */
n1 = 0; /* Position to move remaining elements to */
n2 = 0; /* Position to move remaining elements from */
for (i = 0, api = ATOM_AT_INDEX(mp->atoms, i); i < mp->natoms; i++, api = ATOM_NEXT(api)) {
int j;
Int *cp;
+ cp = AtomConnectData(&api->connect);
for (j = 0; j < api->connect.count; j++) {
- cp = AtomConnectData(&api->connect);
if (cp[j] >= pos)
cp[j]++;
}
+ if (api->anchor != NULL) {
+ cp = AtomConnectData(&api->anchor->connect);
+ for (j = 0; j < api->anchor->connect.count; j++) {
+ if (cp[j] >= pos)
+ cp[j]++;
+ }
+ }
}
for (i = 0; i < mp->nbonds * 2; i++) {
if (mp->bonds[i] >= pos)
if (mp->impropers[i] >= pos)
mp->impropers[i]++;
}
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(mp);
-#endif
}
mp->nframes = -1; /* Should be recalculated later */
MoleculeIncrementModifyCount(mp);
MoleculeCheckSanity(Molecule *mol)
{
const char *fail = "Sanity check failure";
- Int i, j, *ip;
+ Int i, j, *ip, c[4];
Atom *ap;
s_error_count = 0;
for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
ip = AtomConnectData(&ap->connect);
for (j = 0; j < ap->connect.count; j++) {
if (ip[j] < 0 || ip[j] >= mol->natoms)
- s_fprintf(stderr, "%s: atom %d connect %d = %d out of range\n", fail, i, j, ip[j]);
+ s_fprintf(stderr, "%s: atom %d connect[%d] = %d out of range\n", fail, i, j, ip[j]);
if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[j])->connect), i) == 0)
- s_fprintf(stderr, "%s: atom %d connect %d but atom %d has no connect %d\n", fail, i, ip[j], ip[j], i);
+ s_fprintf(stderr, "%s: atom %d has connect %d but atom %d has no connect %d\n", fail, i, ip[j], ip[j], i);
}
}
for (i = 0, ip = mol->bonds; i < mol->nbonds; i++, ip += 2) {
for (i = 0, ip = mol->angles; i < mol->nangles; i++, ip += 3) {
if (ip[0] < 0 || ip[0] >= mol->natoms || ip[1] < 0 || ip[1] >= mol->natoms || ip[2] < 0 || ip[2] >= mol->natoms)
s_fprintf(stderr, "%s: angle %d %d-%d-%d out of range\n", fail, i, ip[0], ip[1], ip[2]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[0])->connect), ip[1]) == 0)
- s_fprintf(stderr, "%s: angle %d %d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[0], ip[1]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[2])->connect), ip[1]) == 0)
- s_fprintf(stderr, "%s: angle %d %d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[2], ip[1]);
+ c[0] = MoleculeAreAtomsConnected(mol, ip[1], ip[0]);
+ if (c[0] == 0)
+ s_fprintf(stderr, "%s: angle %d %d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[1], ip[0]);
+ c[1] = MoleculeAreAtomsConnected(mol, ip[1], ip[2]);
+ if (c[1] == 0)
+ s_fprintf(stderr, "%s: angle %d %d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[1], ip[2]);
+ if (c[0] == 2 && c[1] == 2)
+ s_fprintf(stderr, "%s: angle %d %d-%d-%d but bonds %d-%d and %d-%d are both virtual\n", fail, i, ip[0], ip[1], ip[2], ip[1], ip[0], ip[1], ip[2]);
}
for (i = 0, ip = mol->dihedrals; i < mol->ndihedrals; i++, ip += 4) {
if (ip[0] < 0 || ip[0] >= mol->natoms || ip[1] < 0 || ip[1] >= mol->natoms || ip[2] < 0 || ip[2] >= mol->natoms || ip[3] < 0 || ip[3] >= mol->natoms)
s_fprintf(stderr, "%s: dihedral %d %d-%d-%d%d out of range\n", fail, i, ip[0], ip[1], ip[2], ip[3]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[0])->connect), ip[1]) == 0)
- s_fprintf(stderr, "%s: dihedral %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[0], ip[1]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[1])->connect), ip[2]) == 0)
+ c[0] = MoleculeAreAtomsConnected(mol, ip[1], ip[0]);
+ c[1] = MoleculeAreAtomsConnected(mol, ip[1], ip[2]);
+ c[2] = MoleculeAreAtomsConnected(mol, ip[2], ip[3]);
+ if (c[0] == 0)
+ s_fprintf(stderr, "%s: dihedral %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[1], ip[0]);
+ if (c[1] == 0)
s_fprintf(stderr, "%s: dihedral %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[1], ip[2]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[2])->connect), ip[3]) == 0)
+ if (c[2] == 0)
s_fprintf(stderr, "%s: dihedral %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[2], ip[3]);
}
for (i = 0, ip = mol->impropers; i < mol->nimpropers; i++, ip += 4) {
if (ip[0] < 0 || ip[0] >= mol->natoms || ip[1] < 0 || ip[1] >= mol->natoms || ip[2] < 0 || ip[2] >= mol->natoms || ip[3] < 0 || ip[3] >= mol->natoms)
s_fprintf(stderr, "%s: improper %d %d-%d-%d%d out of range\n", fail, i, ip[0], ip[1], ip[2], ip[3]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[0])->connect), ip[2]) == 0)
- s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[0], ip[2]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[1])->connect), ip[2]) == 0)
- s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[1], ip[2]);
- if (AtomConnectHasEntry(&(ATOM_AT_INDEX(mol->atoms, ip[3])->connect), ip[2]) == 0)
- s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[3], ip[2]);
+ c[0] = MoleculeAreAtomsConnected(mol, ip[2], ip[0]);
+ c[1] = MoleculeAreAtomsConnected(mol, ip[2], ip[1]);
+ c[2] = MoleculeAreAtomsConnected(mol, ip[2], ip[3]);
+ if (c[0] == 0)
+ s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[2], ip[0]);
+ if (c[1] == 0)
+ s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[2], ip[1]);
+ if (c[2] == 0)
+ s_fprintf(stderr, "%s: improper %d %d-%d-%d-%d but atom %d has no connect %d\n", fail, i, ip[0], ip[1], ip[2], ip[3], ip[2], ip[3]);
}
return s_error_count;
}
act = NULL;
__MoleculeLock(dst);
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(dst);
-#endif
+
nsrc = src->natoms;
ndst = dst->natoms;
if (resSeqOffset < 0)
cp[j] = old2new[cp[j] + n1];
if (SYMOP_ALIVE(ap->symop))
ap->symbase = old2new[ap->symbase + n1];
+ if (ap->anchor != NULL) {
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (j = 0; j < ap->anchor->connect.count; j++)
+ cp[j] = old2new[cp[j] + n1];
+ }
}
/* Move the bonds, angles, dihedrals, impropers */
goto panic;
/* Copy the items */
memmove(*items + n1 * nsize, *items_src, sizeof(Int) * nsize * n2);
+ if (i == 0) {
+ /* Copy the bond order info if present */
+ Int nn1 = dst->nbondOrders;
+ if (dst->bondOrders != NULL || src->bondOrders != NULL) {
+ if (AssignArray(&dst->bondOrders, &dst->nbondOrders, sizeof(Double), dst->nbonds - 1, NULL) == NULL)
+ goto panic;
+ memset(dst->bondOrders + nn1, 0, sizeof(Double) * (dst->nbonds - nn1));
+ if (src->bondOrders != NULL)
+ memmove(dst->bondOrders + n1, src->bondOrders, sizeof(Double) * n2);
+ }
+ }
}
/* Renumber */
for (j = 0; j < n1 * nsize; j++)
}
}
-#if PIATOM
- /* Renumber the existing pi-atoms */
- if (dst->npiatoms > 0) {
- for (i = 0; i < dst->npiatoms; i++) {
- PiAtom *pp;
- pp = &dst->piatoms[i];
- cp = AtomConnectData(&pp->connect);
- for (j = 0; j < pp->connect.count; j++) {
- cp[j] = old2new[cp[j]];
- }
- }
- if (dst->npibonds > 0) {
- cp = dst->pibonds;
- for (i = 0; i < dst->npibonds * 4; i++) {
- if (cp[i] < 0)
- continue;
- else if (cp[i] >= ATOMS_MAX_NUMBER)
- continue;
- else {
- cp[i] = old2new[cp[i]];
- }
- }
- }
- }
-
- /* Copy the pi-atoms */
- if (src->npiatoms > 0 && forUndo == 0) {
- int nsrcp, ndstp;
- nsrcp = src->npiatoms;
- ndstp = dst->npiatoms;
- if (AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), nsrcp + ndstp - 1, NULL) == NULL)
- goto panic;
- for (i = 0; i < nsrcp; i++) {
- PiAtom *pp;
- pp = &dst->piatoms[ndstp + i];
- PiAtomDuplicate(pp, &src->piatoms[i]);
- cp = AtomConnectData(&pp->connect);
- for (j = 0; j < pp->connect.count; j++) {
- cp[j] = old2new[ndst + cp[j]];
- }
- if (nactions != NULL) {
- /* This is very inefficient, yet should not cause big problems
- because the number of piatoms in the molecule is usually limited. */
- act = MolActionNew(gMolActionRemoveOnePiAtom, ndstp + i);
- AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
- act = NULL;
- }
- }
- if (src->npibonds > 0) {
- n1 = src->npibonds;
- n2 = dst->npibonds;
- if (AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, n1 + n2 - 1, NULL) == NULL)
- goto panic;
- cp = &dst->pibonds[n2 * 4];
- memmove(cp, src->pibonds, sizeof(Int) * 4 * n1);
- for (i = 0; i < 4 * n1; i++) {
- /* Renumber the pi-atom constructs */
- if (cp[i] < 0)
- continue;
- else if (cp[i] < ATOMS_MAX_NUMBER)
- cp[i] = old2new[ndst + cp[i]]; /* Ordinary atoms */
- else
- cp[i] += ndstp; /* pi-atoms */
- }
- if (nactions != NULL) {
- ig = IntGroupNewWithPoints(n2, n1, -1);
- act = MolActionNew(gMolActionRemovePiBonds, ig);
- AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
- act = NULL;
- IntGroupRelease(ig);
- }
- }
- }
-#endif /* PIATOM */
-
- MoleculeCleanUpResidueTable(dst);
+ MoleculeCleanUpResidueTable(dst);
free(new2old);
dst->nframes = -1; /* Should be recalculated later */
return 1; /* Not reached */
}
+/* Unmerge the molecule. If necessary, the undo actions are stored in nactions/actions array.
+ (The nactions/actions array must be initialized by the caller) */
static int
sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqOffset, int moveFlag, Int *nactions, MolAction ***actions, Int forUndo)
{
return 1; /* Bad parameter */
__MoleculeLock(src);
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(src);
-#endif
- if (nactions != NULL)
- *nactions = 0;
- if (actions != NULL)
- *actions = NULL;
act = NULL;
nsrc = src->natoms;
}
}
} else dst_par_g = remove_par_g = NULL;
-
+
+ /* Pi anchors should be modified if the anchor and its component atoms become separated between
+ src anc dst */
+ if (moveFlag) {
+ Int ibufsize, *ibuf, flag_i, flag_j;
+ ibufsize = 8;
+ ibuf = (Int *)malloc(sizeof(Int) * ibufsize);
+ for (i = 0, ap = src->atoms; i < src->natoms; i++, ap = ATOM_NEXT(ap)) {
+ if (ap->anchor == NULL)
+ continue;
+ flag_i = (old2new[i] < nsrcnew);
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (j = n1 = 0; j < ap->anchor->connect.count; j++) {
+ flag_j = (old2new[cp[j]] < nsrcnew);
+ if (flag_i == flag_j) {
+ if (n1 >= ibufsize) {
+ ibufsize += 8;
+ ibuf = (Int *)realloc(ibuf, sizeof(Int) * ibufsize);
+ }
+ ibuf[n1++] = cp[j];
+ }
+ }
+ if (n1 < j) {
+ /* Need to modify the pi anchor list */
+ if (n1 <= 1)
+ n1 = 0;
+ MolActionCreateAndPerform(src, SCRIPT_ACTION("isI"), "set_atom_attr", i, "anchor_list", n1, ibuf);
+ }
+ }
+ }
+
/* Make a new molecule */
if (dstp != NULL) {
dst = MoleculeNew();
dst_ap = NULL;
}
- /* Renumber the atom indices in connect[] */
+ /* Renumber the atom indices in connect[] (src) */
if (moveFlag) {
for (i = 0, ap = src->atoms; i < src->natoms; i++, ap = ATOM_NEXT(ap)) {
cp = AtomConnectData(&ap->connect);
cp[n1++] = n2;
}
AtomConnectResize(&ap->connect, n1);
+ if (ap->anchor != NULL) {
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (j = n1 = 0; j < ap->anchor->connect.count; j++) {
+ n2 = old2new[cp[j]];
+ if (n2 < nsrcnew)
+ cp[n1++] = n2;
+ }
+ if (n1 != ap->anchor->connect.count) {
+ /* This should not happen!! */
+ AtomConnectResize(&ap->anchor->connect, n1);
+ fprintf(stderr, "Internal error in sMoleculeUnmergeSub (line %d)\n", __LINE__);
+ if (n1 == 0) {
+ free(ap->anchor->coeffs);
+ free(ap->anchor);
+ ap->anchor = NULL;
+ }
+ }
+ }
}
}
- /* Renumber the atom indices in connect[] and the residue indices */
+ /* Renumber the atom indices in connect[] (dst) */
if (dst != NULL) {
for (i = 0, ap = dst->atoms; i < dst->natoms; i++, ap = ATOM_NEXT(ap)) {
if (ap->resSeq != 0 && ap->resSeq - resSeqOffset >= 0)
cp[n1++] = n2;
}
AtomConnectResize(&ap->connect, n1);
+ if (ap->anchor != NULL) {
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (j = n1 = 0; j < ap->anchor->connect.count; j++) {
+ n2 = old2new[cp[j]] - nsrcnew;
+ if (n2 >= 0)
+ cp[n1++] = n2;
+ }
+ if (n1 != ap->anchor->connect.count) {
+ /* This can happen, and the anchor info is silently modified */
+ if (n1 <= 1) {
+ AtomConnectResize(&ap->anchor->connect, 0);
+ free(ap->anchor->coeffs);
+ free(ap->anchor);
+ ap->anchor = NULL;
+ } else {
+ Double d;
+ AtomConnectResize(&ap->anchor->connect, n1);
+ d = 0.0;
+ for (j = 0; j < n1; j++)
+ d += ap->anchor->coeffs[j];
+ for (j = 0; j < n1; j++)
+ ap->anchor->coeffs[j] /= d;
+ MoleculeCalculatePiAnchorPosition(dst, i);
+ }
+ }
+ }
}
}
goto panic;
if (sCopyElementsFromArrayAtPositions(*items, *nitems, *items_dst, sizeof(Int) * nsize, move_g) != 0)
goto panic;
+ if (i == 0 && src->bondOrders != NULL) {
+ if (AssignArray(&dst->bondOrders, &dst->nbondOrders, sizeof(Double), n3 - 1, NULL) == NULL)
+ goto panic;
+ if (sCopyElementsFromArrayAtPositions(src->bondOrders, src->nbondOrders, dst->bondOrders, sizeof(Double), move_g) != 0)
+ goto panic;
+ }
}
/* Remove from src */
if (moveFlag && forUndo == 0) {
if (nactions != NULL) {
Int k, *ip;
+ Double *dp;
ip = (Int *)malloc(sizeof(Int) * nsize * n2);
for (j = 0; (k = IntGroupGetNthPoint(del_g, j)) >= 0; j++)
memmove(ip + j * nsize, *items + k * nsize, sizeof(Int) * nsize);
+ if (i == 0 && src->bondOrders != NULL) {
+ dp = (Double *)malloc(sizeof(Double) * n2);
+ for (j = 0; (k = IntGroupGetNthPoint(del_g, j)) >= 0; j++)
+ dp[j] = src->bondOrders[k];
+ } else dp = NULL;
switch (i) {
case 0:
act = MolActionNew(gMolActionAddBondsForUndo, n2 * nsize, ip, del_g); break;
act = NULL;
}
free(ip);
+ if (dp != NULL) {
+ act = MolActionNew(gMolActionAssignBondOrders, n2, dp, del_g);
+ AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
+ act = NULL;
+ free(dp);
+ }
}
if (sRemoveElementsFromArrayAtPositions(*items, *nitems, NULL, sizeof(Int) * nsize, del_g) != 0)
goto panic;
}
}
-#if PIATOM
- /* Copy the pi-atoms */
- if (src->npiatoms > 0) {
- PiAtom *pp1, *pp2;
- Int *patoms_old2new;
- n1 = src->npiatoms;
- patoms_old2new = (Int *)calloc(sizeof(Int), src->npiatoms);
- if (patoms_old2new == NULL)
- goto panic;
- for (i = 0, pp1 = src->piatoms; i < src->npiatoms; i++, pp1++) {
- /* Is this entry to be copied to dst? */
- cp = AtomConnectData(&pp1->connect);
- for (j = pp1->connect.count - 1; j >= 0; j--) {
- if (old2new[cp[j]] < nsrcnew)
- break;
- }
- if (j < 0) {
- /* Copy this entry */
- patoms_old2new[i] = dst->npiatoms;
- pp2 = AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), dst->npiatoms, NULL);
- PiAtomDuplicate(pp2, pp1);
- cp = AtomConnectData(&pp2->connect);
- for (j = 0; j < pp2->connect.count; j++)
- cp[j] = old2new[cp[j]] - nsrcnew;
- } else {
- patoms_old2new[i] = -1;
- }
- }
- /* Copy the piatom constructs to dst */
- for (i = 0; i < src->npibonds; i++) {
- for (j = 0; j < 4; j++) {
- n2 = src->pibonds[i * 4 + j];
- if (n2 >= ATOMS_MAX_NUMBER) {
- if (patoms_old2new[n2 - ATOMS_MAX_NUMBER] < 0)
- break;
- } else if (n2 >= 0) {
- if (old2new[n2] < nsrcnew)
- break;
- }
- }
- if (j >= 4) {
- /* Copy this entry */
- cp = (Int *)AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, dst->npibonds, &src->pibonds[i * 4]);
- for (j = 0; j < 4; j++) {
- if (cp[j] >= ATOMS_MAX_NUMBER) {
- cp[j] = patoms_old2new[cp[j] - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
- } else if (cp[j] >= 0) {
- cp[j] = old2new[cp[j]] - nsrcnew;
- }
- }
- }
- }
- if (moveFlag) {
- Int npibonds_to_go, *pibonds_to_go;
- IntGroup *pibonds_group;
-
- /* Remove the piatom entries containing non-remaining atoms. Note: the piatom
- entries that do not remain in src and not copied to dst will disappear. */
- n2 = 0;
- for (i = 0; i < src->npiatoms; i++) {
- /* Is this entry to be removed? */
- pp1 = &src->piatoms[i];
- cp = AtomConnectData(&pp1->connect);
- for (j = pp1->connect.count - 1; j >= 0; j--) {
- if (old2new[cp[j]] >= nsrcnew)
- break;
- }
- patoms_old2new[i] = (j < 0 ? n2++ : -1);
- }
-
- /* Remove pibonds first (necessary for undo) */
- npibonds_to_go = 0;
- pibonds_to_go = NULL;
- pibonds_group = NULL;
- for (i = src->npibonds - 1; i >= 0; i--) {
- for (j = 0; j < 4; j++) {
- n3 = src->pibonds[i * 4 + j];
- if (n3 >= ATOMS_MAX_NUMBER) {
- if (patoms_old2new[n3 - ATOMS_MAX_NUMBER] < 0)
- break;
- } else if (n3 >= 0) {
- if (old2new[n3] >= nsrcnew)
- break;
- }
- }
- if (j < 4) {
- /* Remove */
- if (nactions != NULL) {
- /* Since we are scanning pibonds[] from the end, the new entry should be inserted
- at the top of the pibonds_to_go[] array. */
- InsertArray(&pibonds_to_go, &npibonds_to_go, sizeof(Int) * 4, 0, 1, src->pibonds + i * 4);
- if (pibonds_group == NULL)
- pibonds_group = IntGroupNew();
- IntGroupAdd(pibonds_group, i, 1);
- }
- DeleteArray(&src->pibonds, &src->npibonds, sizeof(Int) * 4, i, 1, NULL);
- } else {
- /* Renumber */
- cp = &src->pibonds[i * 4];
- for (j = 0; j < 4; j++) {
- n3 = cp[j];
- if (n3 >= ATOMS_MAX_NUMBER)
- cp[j] = patoms_old2new[n3 - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
- else if (n3 >= 0)
- cp[j] = old2new[n3];
- }
- }
- }
- if (nactions != NULL && pibonds_to_go != NULL) {
- act = MolActionNew(gMolActionInsertPiBonds, pibonds_group, npibonds_to_go * 4, pibonds_to_go);
- AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
- act = NULL;
- free(pibonds_to_go);
- }
- if (pibonds_group != NULL)
- IntGroupRelease(pibonds_group);
-
- for (i = src->npiatoms - 1; i >= 0; i--) {
- pp1 = &src->piatoms[i];
- if (patoms_old2new[i] < 0) {
- /* Remove the entries */
- /* (If forUndo is true, these entries should already have been removed.) */
- if (nactions != NULL) {
- act = MolActionNew(gMolActionInsertOnePiAtom, i, 4, pp1->aname, pp1->type, pp1->connect.count, AtomConnectData(&pp1->connect), pp1->ncoeffs, pp1->coeffs);
- AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
- act = NULL;
- }
- PiAtomClean(pp1);
- DeleteArray(&src->piatoms, &src->npiatoms, sizeof(PiAtom), i, 1, NULL);
- } else {
- /* Renumber */
- cp = AtomConnectData(&pp1->connect);
- for (j = 0; j < pp1->connect.count; j++) {
- cp[j] = old2new[cp[j]];
- }
- }
- }
- }
-
- free(patoms_old2new);
- }
-#endif /* PIATOM */
-
/* Clean up */
IntGroupRelease(remain_g);
MoleculeCleanUpResidueTable(src);
int
MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, Int autoGenerate)
{
- int i, j, n1, n2;
- Atom *ap;
- Int *cp;
-
+ Int nangles, ndihedrals;
+ Int *angles, *dihedrals;
+ Int i, j, k, kk, n1, n2, cn1, cn2;
+ Int *cp1, *cp2;
+ Int temp[4];
+ Atom *ap1, *ap2, *ap3;
+
if (mp == NULL || bonds == NULL || nbonds <= 0)
return 0;
if (mp->noModifyTopology)
__MoleculeUnlock(mp);
return -4; /* Out of memory */
}
+ if (mp->bondOrders != NULL) {
+ /* Expand the bond order info (all new entries are zero) */
+ Double *dp = (Double *)calloc(sizeof(Double), nbonds);
+ if (dp == NULL)
+ return -4;
+ if (AssignArray(&(mp->bondOrders), &(mp->nbondOrders), sizeof(Double), n1 + nbonds - 1, NULL) == NULL
+ || sInsertElementsToArrayAtPositions(mp->bondOrders, n1, dp, nbonds, sizeof(Double), where) != 0) {
+ __MoleculeUnlock(mp);
+ free(dp);
+ return -4;
+ }
+ free(dp);
+ }
+
+ angles = dihedrals = NULL;
+ nangles = ndihedrals = 0;
- /* Add connects[] */
+ /* Add connects[], and angles/dihedrals (if autoGenerate is true) */
for (i = 0; i < nbonds; i++) {
+
+ /* One entry at time */
+ /* (Otherwise, duplicate entries of angles and dihedrals result) */
n1 = bonds[i * 2];
n2 = bonds[i * 2 + 1];
- ap = ATOM_AT_INDEX(mp->atoms, n1);
- cp = AtomConnectData(&ap->connect);
- AtomConnectInsertEntry(&ap->connect, -1, n2);
- ap = ATOM_AT_INDEX(mp->atoms, n2);
- cp = AtomConnectData(&ap->connect);
- AtomConnectInsertEntry(&ap->connect, -1, n1);
- }
+
+ ap1 = ATOM_AT_INDEX(mp->atoms, n1);
+ AtomConnectInsertEntry(&ap1->connect, -1, n2);
+ ap2 = ATOM_AT_INDEX(mp->atoms, n2);
+ AtomConnectInsertEntry(&ap2->connect, -1, n1);
- /* Add angles, dihedrals, impropers */
- if (autoGenerate) {
- Int nangles, ndihedrals;
- Int *angles, *dihedrals;
- Int k, kk;
- Int *cp1, *cp2;
- Int temp[4];
- Atom *ap1, *ap2, *ap3;
-
- angles = dihedrals = NULL;
- nangles = ndihedrals = 0;
-
- for (i = 0; i < nbonds; i++) {
+ /* Add angles and dihedrals */
+ if (autoGenerate) {
AtomConnect *ac1, *ac2;
- n1 = bonds[i * 2];
- n2 = bonds[i * 2 + 1];
- ap1 = ATOM_AT_INDEX(mp->atoms, n1);
- ap2 = ATOM_AT_INDEX(mp->atoms, n2);
if (ap1->anchor == NULL || ap2->anchor == NULL) {
/* N1-N2-{XY} or N2-N1-{XY} angles (X: connected atom, Y: constitute atom of pi-anchor) */
for (j = 0; j < 4; j++) {
case 3: if (ap1->anchor == NULL) continue; else ac1 = &ap1->anchor->connect; break; /* N2-N1-Y */
}
cp1 = AtomConnectData(ac1);
- for (k = 0; k < ac1->count; k++) {
+ cn1 = ac1->count;
+ for (k = 0; k < cn1; k++) {
temp[2] = cp1[k];
if (temp[2] == temp[0])
continue;
+ ap3 = ATOM_AT_INDEX(mp->atoms, temp[2]);
+ if (ap3->anchor != NULL) {
+ /* Avoid X-anchor-anchor angle (anchor-X-anchor is allowed) */
+ if ((j < 2 && ap2->anchor != NULL) || (j >= 2 && ap1->anchor != NULL))
+ continue;
+ }
if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
goto panic;
- /* Dihedrals N1-N2-X-X or N2-N1-X-X */
+ /* Dihedrals N1-N2-X-{XY} or N2-N1-X-{XY} */
if (j == 1 || j == 3)
continue;
- ap3 = ATOM_AT_INDEX(mp->atoms, temp[2]);
cp2 = AtomConnectData(&ap3->connect);
for (kk = 0; kk < ap3->connect.count; kk++) {
temp[3] = cp2[kk];
if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
goto panic;
}
+ if (ap3->anchor != NULL) {
+ /* N1-N2-X-Y or N2-N1-X-Y */
+ /* for Y, only the first constitute atom is considered */
+ cp2 = AtomConnectData(&ap3->anchor->connect);
+ temp[3] = cp2[0];
+ if (temp[3] == temp[0] || temp[3] == temp[1])
+ continue;
+ if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+ goto panic;
+ }
}
}
}
- /* X-N1-N2-X angles */
- if (ap1->anchor == NULL)
+ /* X-N1-N2-X dihedrals */
+ /* Y-N1-N2-anchor is allowed, but the force may be zero if the angle N1-N2-anchor is */
+ /* close to 180 deg (e.g. in ferrocene, C-anchor-Fe-anchor dihedral should be k=0) */
+ if (ap1->anchor == NULL) {
ac1 = &ap1->connect;
- else ac1 = &ap1->anchor->connect;
- if (ap2->anchor == NULL)
+ cn1 = ac1->count;
+ } else {
+ ac1 = &ap1->anchor->connect;
+ cn1 = 1; /* Only the first constitute atom of pi-anchor is considered */
+ }
+ if (ap2->anchor == NULL) {
ac2 = &ap2->connect;
- else ac2 = &ap2->anchor->connect;
+ cn2 = ac2->count;
+ } else {
+ ac2 = &ap2->anchor->connect;
+ cn2 = 1; /* Only the first constitute atom of pi-anchor is considered */
+ }
temp[1] = n1;
temp[2] = n2;
cp1 = AtomConnectData(ac1);
cp2 = AtomConnectData(ac2);
- for (j = 0; j < ac1->count; j++) {
+ for (j = 0; j < cn1; j++) {
temp[0] = cp1[j];
if (temp[0] == temp[2])
continue;
- if (ATOM_AT_INDEX(mp->atoms, temp[0])->anchor != NULL)
- continue;
- for (k = 0; k < ac2->count; k++) {
+ for (k = 0; k < cn2; k++) {
temp[3] = cp2[k];
if (temp[3] == temp[0] || temp[3] == temp[1])
continue;
- if (ATOM_AT_INDEX(mp->atoms, temp[3])->anchor != NULL)
- continue;
if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
goto panic;
}
}
}
- temp[0] = kInvalidIndex; /* For termination */
- if (angles != NULL) {
- if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
- goto panic;
- MoleculeAddAngles(mp, angles, NULL);
- free(angles);
- }
- if (dihedrals != NULL) {
- if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
- goto panic;
- MoleculeAddDihedrals(mp, dihedrals, NULL);
- free(dihedrals);
- }
}
+ if (angles != NULL) {
+ temp[0] = kInvalidIndex;
+ if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
+ goto panic;
+ MoleculeAddAngles(mp, angles, NULL);
+ free(angles);
+ }
+ if (dihedrals != NULL) {
+ temp[0] = kInvalidIndex;
+ if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+ goto panic;
+ MoleculeAddDihedrals(mp, dihedrals, NULL);
+ free(dihedrals);
+ }
+
MoleculeIncrementModifyCount(mp);
mp->needsMDRebuild = 1;
__MoleculeUnlock(mp);
int
MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved, IntGroup **outRemovedPos)
{
- Int i, j, n1, n2;
- Int *ip, na, nd, ni;
+ Int i, j, n1, n2, nw;
+ Int *ip, *jp, na, nd, ni;
IntGroup *ag, *dg, *ig;
Atom *ap;
IntGroupIterator iter;
}
}
}
- IntGroupIteratorReset(&iter);
/* Remove bonds, angles, dihedrals, impropers */
- ag = dg = ig = NULL;
+ ag = IntGroupNew();
+ dg = ig = NULL;
na = nd = ni = 0;
+
+ nw = IntGroupGetCount(where);
+ jp = (Int *)malloc(sizeof(Int) * nw * 2);
+ j = 0;
+ IntGroupIteratorReset(&iter);
while ((i = IntGroupIteratorNext(&iter)) >= 0) {
- n1 = mp->bonds[i * 2];
- n2 = mp->bonds[i * 2 + 1];
- for (j = 0; j < mp->nangles; j++) {
- ip = mp->angles + j * 3;
+ jp[j++] = mp->bonds[i * 2];
+ jp[j++] = mp->bonds[i * 2 + 1];
+ }
+ IntGroupIteratorRelease(&iter);
+
+ for (i = 0, ip = mp->angles; i < mp->nangles; i++, ip += 3) {
+ for (j = 0; j < nw; j++) {
+ n1 = jp[j * 2];
+ n2 = jp[j * 2 + 1];
if ((ip[0] == n1 && ip[1] == n2)
- || (ip[1] == n1 && ip[0] == n2)
- || (ip[1] == n1 && ip[2] == n2)
- || (ip[2] == n1 && ip[1] == n2)) {
- if (ag == NULL)
- ag = IntGroupNew();
- if (IntGroupAdd(ag, j, 1) != 0)
+ || (ip[1] == n1 && ip[0] == n2)
+ || (ip[1] == n1 && ip[2] == n2)
+ || (ip[2] == n1 && ip[1] == n2)) {
+ if (IntGroupAdd(ag, i, 1) != 0)
goto panic;
na++;
+ break;
}
}
- for (j = 0; j < mp->ndihedrals; j++) {
- ip = mp->dihedrals + j * 4;
+ }
+ for (i = 0, ip = mp->dihedrals; i < mp->ndihedrals; i++, ip += 4) {
+ for (j = 0; j < nw; j++) {
+ n1 = jp[j * 2];
+ n2 = jp[j * 2 + 1];
if ((ip[0] == n1 && ip[1] == n2)
|| (ip[1] == n1 && ip[0] == n2)
|| (ip[1] == n1 && ip[2] == n2)
|| (ip[3] == n1 && ip[2] == n2)) {
if (dg == NULL)
dg = IntGroupNew();
- if (IntGroupAdd(dg, j, 1) != 0)
+ if (IntGroupAdd(dg, i, 1) != 0)
goto panic;
nd++;
+ break;
}
}
- for (j = 0; j < mp->nimpropers; j++) {
- ip = mp->impropers + j * 4;
+ }
+ for (i = 0, ip = mp->impropers; i < mp->nimpropers; i++, ip += 4) {
+ for (j = 0; j < nw; j++) {
+ n1 = jp[j * 2];
+ n2 = jp[j * 2 + 1];
if ((ip[0] == n1 && ip[2] == n2)
|| (ip[1] == n1 && ip[2] == n2)
|| (ip[3] == n1 && ip[2] == n2)
|| (ip[3] == n2 && ip[2] == n1)) {
if (ig == NULL)
ig = IntGroupNew();
- if (IntGroupAdd(ig, j, 1) != 0)
+ if (IntGroupAdd(ig, i, 1) != 0)
goto panic;
ni++;
+ break;
}
}
}
- IntGroupIteratorRelease(&iter);
-
+ free(jp);
+
if (sRemoveElementsFromArrayAtPositions(mp->bonds, mp->nbonds, bonds, sizeof(Int) * 2, where) != 0)
goto panic;
mp->nbonds -= IntGroupGetCount(where);
free(mp->bonds);
mp->bonds = NULL;
}
+ if (mp->bondOrders != NULL) {
+ if (sRemoveElementsFromArrayAtPositions(mp->bondOrders, mp->nbondOrders, NULL, sizeof(Double), where) != 0)
+ goto panic;
+ mp->nbondOrders -= IntGroupGetCount(where);
+ if (mp->nbondOrders == 0) {
+ free(mp->bondOrders);
+ mp->bondOrders = NULL;
+ }
+ }
if (na == 0 && nd == 0 && ni == 0)
ip = NULL;
else
IntGroupRelease(ig);
}
+ if (IntGroupGetCount(ag) == 0) {
+ IntGroupRelease(ag);
+ ag = NULL;
+ }
+
*outRemoved = ip;
*outRemovedPos = ag;
}
int
+MoleculeAssignBondOrders(Molecule *mp, const Double *orders, IntGroup *where)
+{
+ Int i, j;
+ IntGroupIterator iter;
+ if (mp == NULL || orders == NULL || mp->nbonds == 0)
+ return 0;
+ if (mp->noModifyTopology)
+ return -4; /* Prohibited operation */
+ if (mp->bondOrders == NULL) {
+ AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbonds - 1, NULL);
+ memset(mp->bondOrders, 0, sizeof(Double) * mp->nbondOrders);
+ }
+ IntGroupIteratorInit(where, &iter);
+ j = 0;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ if (i >= mp->nbondOrders)
+ break;
+ mp->bondOrders[i] = orders[j++];
+ }
+ IntGroupIteratorRelease(&iter);
+ return 0;
+}
+
+int
+MoleculeGetBondOrders(Molecule *mp, Double *outOrders, IntGroup *where)
+{
+ Int i, j;
+ IntGroupIterator iter;
+ if (mp == NULL || mp->nbonds == 0)
+ return 0;
+ if (mp->bondOrders == NULL) {
+ /* Returns all zero */
+ i = IntGroupGetCount(where);
+ for (j = 0; j < i; j++)
+ outOrders[j] = 0.0;
+ } else {
+ IntGroupIteratorInit(where, &iter);
+ j = 0;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ if (i < mp->nbondOrders)
+ outOrders[j] = mp->bondOrders[i];
+ else outOrders[j] = 0.0;
+ j++;
+ }
+ }
+ return 0;
+}
+
+int
MoleculeAddAngles(Molecule *mp, const Int *angles, IntGroup *where)
{
int n1, nc;
__MoleculeLock(mp);
if (sRemoveElementsFromArrayAtPositions(mp->angles, mp->nangles, angles, sizeof(Int) * 3, where) != 0) {
__MoleculeUnlock(mp);
- Panic("Low memory while adding angles");
+ Panic("Bad argument while deleting angles");
}
mp->nangles -= (nc = IntGroupGetCount(where));
if (mp->nangles == 0) {
__MoleculeLock(mp);
if (sRemoveElementsFromArrayAtPositions(mp->dihedrals, mp->ndihedrals, dihedrals, sizeof(Int) * 4, where) != 0) {
__MoleculeUnlock(mp);
- Panic("Low memory while adding dihedrals");
+ Panic("Internal error: bad argument while deleting dihedrals");
}
mp->ndihedrals -= (nc = IntGroupGetCount(where));
if (mp->ndihedrals == 0) {
__MoleculeLock(mp);
if (sRemoveElementsFromArrayAtPositions(mp->impropers, mp->nimpropers, impropers, sizeof(Int) * 4, where) != 0) {
__MoleculeUnlock(mp);
- Panic("Low memory while adding impropers");
+ Panic("Internal error: bad argument while deleting impropers");
}
mp->nimpropers -= (nc = IntGroupGetCount(where));
if (mp->impropers == NULL) {
Atom *ap;
/* If LSB of resSeqs is 1, then a constant value is used for all specified atoms */
- if (((int)resSeqs & 1) == 0) {
+ if (((uintptr_t)resSeqs & 1) == 0) {
withArray = 1;
resSeq = 0;
} else {
withArray = 0;
- resSeq = ((int)resSeqs - 1) / 2;
+ resSeq = ((uintptr_t)resSeqs - 1) / 2;
}
IntGroupIteratorInit(group, &iter);
int
MoleculeChangeResidueNumber(Molecule *mp, IntGroup *group, int resSeq)
{
- return MoleculeChangeResidueNumberWithArray(mp, group, (Int *)(resSeq * 2 + 1));
+ return MoleculeChangeResidueNumberWithArray(mp, group, (Int *)(intptr_t)(resSeq * 2 + 1));
}
/* Offset the residue numbers by a certain amount. The argument nresidues, if non-negative,
free(newAtoms);
free(old2new);
free(apArray);
-
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(mp);
-#endif
}
/* Renumber atoms */
for (i = 0; i < mp->nimpropers * 4; i++) {
mp->impropers[i] = old2new[mp->impropers[i]];
}
+ /* Renumber the connection table and pi anchor table */
for (i = 0; i < mp->natoms; i++) {
Atom *ap = ATOM_AT_INDEX(saveAtoms, i);
Int *ip = AtomConnectData(&ap->connect);
for (j = 0; j < ap->connect.count; j++, ip++)
*ip = old2new[*ip];
+ if (ap->anchor != NULL) {
+ ip = AtomConnectData(&ap->anchor->connect);
+ for (j = 0; j < ap->anchor->connect.count; j++, ip++)
+ *ip = old2new[*ip];
+ }
}
+
if (mp->par != NULL) {
/* Renumber the parameters */
int n;
}
}
-#if PIATOM
- if (mp->npiatoms) {
- /* Renumber the pi-atoms */
- for (i = 0; i < mp->npiatoms; i++) {
- PiAtom *pp = &mp->piatoms[i];
- Int *cp = AtomConnectData(&pp->connect);
- for (j = 0; j < pp->connect.count; j++)
- cp[j] = old2new[cp[j]];
- }
- }
-
- if (mp->npibonds) {
- /* Renumber the pi-atom constructs */
- for (i = 0; i < mp->npibonds * 4; i++) {
- j = mp->pibonds[i];
- if (j >= 0 && j < mp->natoms)
- mp->pibonds[i] = old2new[j];
- }
- }
-#endif /* PIATOM */
-
/* Renumber the atoms */
for (i = 0; i < mp->natoms; i++)
memmove(ATOM_AT_INDEX(mp->atoms, old2new[i]), ATOM_AT_INDEX(saveAtoms, i), gSizeOfAtomRecord);
retval = 0;
-#if PIATOM
- MoleculeInvalidatePiConnectionTable(mp);
-#endif
MoleculeIncrementModifyCount(mp);
mp->needsMDRebuild = 1;
return -1; /* Non-regular transform */
/* Calculate the reciprocal cell parameters */
- cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[1] * cp->rtr[1] + cp->rtr[2] * cp->rtr[2]);
- cp->rcell[1] = sqrt(cp->rtr[3] * cp->rtr[3] + cp->rtr[4] * cp->rtr[4] + cp->rtr[5] * cp->rtr[5]);
- cp->rcell[2] = sqrt(cp->rtr[6] * cp->rtr[6] + cp->rtr[7] * cp->rtr[7] + cp->rtr[8] * cp->rtr[8]);
- cp->rcell[3] = acos((cp->rtr[3] * cp->rtr[6] + cp->rtr[4] * cp->rtr[7] + cp->rtr[5] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
- cp->rcell[4] = acos((cp->rtr[6] * cp->rtr[0] + cp->rtr[7] * cp->rtr[1] + cp->rtr[8] * cp->rtr[2]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
- cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[3] + cp->rtr[1] * cp->rtr[4] + cp->rtr[2] * cp->rtr[5]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
+ cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[3] * cp->rtr[3] + cp->rtr[6] * cp->rtr[6]);
+ cp->rcell[1] = sqrt(cp->rtr[1] * cp->rtr[1] + cp->rtr[4] * cp->rtr[4] + cp->rtr[7] * cp->rtr[7]);
+ cp->rcell[2] = sqrt(cp->rtr[2] * cp->rtr[2] + cp->rtr[5] * cp->rtr[5] + cp->rtr[8] * cp->rtr[8]);
+ cp->rcell[3] = acos((cp->rtr[1] * cp->rtr[2] + cp->rtr[4] * cp->rtr[5] + cp->rtr[7] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
+ cp->rcell[4] = acos((cp->rtr[2] * cp->rtr[0] + cp->rtr[5] * cp->rtr[3] + cp->rtr[8] * cp->rtr[6]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
+ cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[1] + cp->rtr[3] * cp->rtr[4] + cp->rtr[6] * cp->rtr[7]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
if (calc_abc) {
/* Calculate a, b, c, alpha, beta, gamma */
cp->cell[4] = acos((cp->tr[6] * cp->tr[0] + cp->tr[7] * cp->tr[1] + cp->tr[8] * cp->tr[2]) / (cp->cell[2] * cp->cell[0])) * kRad2Deg;
cp->cell[5] = acos((cp->tr[0] * cp->tr[3] + cp->tr[1] * cp->tr[4] + cp->tr[2] * cp->tr[5]) / (cp->cell[0] * cp->cell[1])) * kRad2Deg;
}
-
return 0;
}
}
MatrixSymDiagonalize(m1, val, axis);
for (u = 0; u < 3; u++) {
+ anp->eigval[u] = val[u];
if (val[u] < 0) {
fprintf(stderr, "Non-positive definite thermal parameters for atom %.4s\n", mp->atoms[n1].aname);
val[u] = 0.001;
idx2 = cp[i];
if (IntGroupLookup(result, idx2, NULL))
continue;
+ if (ap->anchor != NULL && ATOM_AT_INDEX(mp->atoms, idx2)->anchor != NULL)
+ continue; /* bond between two pi_anchors is ignored */
sMoleculeFragmentSub(mp, idx2, result, exatoms);
}
-#if PIATOM
- if (mp->piconnects != NULL) {
- for (i = mp->piconnects[idx]; i < mp->piconnects[idx + 1]; i++) {
- idx2 = mp->piconnects[i];
+ if (ap->anchor != NULL) {
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (i = 0; i < ap->anchor->connect.count; i++) {
+ idx2 = cp[i];
if (IntGroupLookup(result, idx2, NULL))
continue;
sMoleculeFragmentSub(mp, idx2, result, exatoms);
}
}
-#endif
}
/* 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();
-#if PIATOM
- MoleculeValidatePiConnectionTable(mp);
-#endif
sMoleculeFragmentSub(mp, n1, result, exatoms);
return result;
}
for (i = 0; i < argc; i++)
IntGroupAdd(exatoms, argv[i], 1);
result = IntGroupNew();
-#if PIATOM
- MoleculeValidatePiConnectionTable(mp);
-#endif
sMoleculeFragmentSub(mp, n1, result, exatoms);
IntGroupRelease(exatoms);
return result;
return NULL;
IntGroupIteratorInit(inatoms, &iter);
result = IntGroupNew();
-#if PIATOM
- MoleculeValidatePiConnectionTable(mp);
-#endif
while ((i = IntGroupIteratorNext(&iter)) >= 0) {
sMoleculeFragmentSub(mp, i, result, exatoms);
}
Atom *ap;
if (mp == NULL || mp->natoms == 0 || group == NULL)
return 0; /* Invalid arguments */
-#if PIATOM
- MoleculeValidatePiConnectionTable(mp);
-#endif
bond_count = 0;
for (i = 0; (i1 = IntGroupGetStartPoint(group, i)) >= 0; i++) {
i2 = IntGroupGetEndPoint(group, i);
ap = ATOM_AT_INDEX(mp->atoms, j);
cp = AtomConnectData(&ap->connect);
for (k = 0; k < ap->connect.count; k++) {
+ if (ap->anchor != NULL && ATOM_AT_INDEX(mp->atoms, cp[k])->anchor != NULL)
+ continue; /* Ignore bond between two pi_anchors */
if (IntGroupLookup(group, cp[k], NULL) == 0) {
bond_count++;
nval1 = j;
return 0; /* Too many bonds */
}
}
-#if PIATOM
- 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) {
+ if (ap->anchor != NULL) {
+ cp = AtomConnectData(&ap->anchor->connect);
+ for (k = 0; k < ap->anchor->connect.count; k++) {
+ if (IntGroupLookup(group, cp[k], NULL) == 0) {
bond_count++;
nval1 = j;
- nval2 = k2;
+ nval2 = cp[k];
if (bond_count > 1)
return 0; /* Too many bonds */
}
- }
+ }
}
-#endif
}
}
if (bond_count == 1) {
int i, j, count, n_new, n_old, natoms, exframes, last_inserted;
Vector *tempv, *vp;
Atom *ap;
+ MolProp *prp;
+ Double *dp;
+
if (mp == NULL || (natoms = mp->natoms) == 0 || (count = IntGroupGetCount(group)) <= 0)
return -1;
/* Expand ap->frames for all atoms */
for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
- if (ap->frames == NULL)
- vp = (Vector *)calloc(sizeof(Vector), n_new);
- else
- vp = (Vector *)realloc(ap->frames, sizeof(Vector) * n_new);
- if (vp == NULL) {
+ Int n = ap->nframes;
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), n_new - 1, NULL);
+ if (ap->frames == NULL) {
__MoleculeUnlock(mp);
return -1;
}
- for (j = ap->nframes; j < n_new; j++)
- vp[j] = ap->r;
- ap->frames = vp;
+ for (j = n; j < n_new; j++)
+ ap->frames[j] = ap->r;
}
if (mp->cell != NULL) {
j = mp->nframe_cells;
}
}
+ /* Expand propvals for all properties */
+ for (i = 0, prp = mp->molprops; i < mp->nmolprops; i++, prp++) {
+ dp = (Double *)realloc(prp->propvals, sizeof(Double) * n_new);
+ if (dp == NULL) {
+ __MoleculeUnlock(mp);
+ return -1;
+ }
+ for (j = n_old; j < n_new; j++)
+ dp[j] = 0.0;
+ prp->propvals = dp;
+ }
+
/* group = [n0..n1-1, n2..n3-1, ...] */
/* s = t = 0, */
/* tempv[0..n0-1] <- ap[0..n0-1], s += n0,
...
tempv[nl..n_new-1] <- ap[s..s+(n_new-nl-1)], s += n_new-nl
At last, s will become n_old and t will become count. */
- for (i = 0, ap = mp->atoms; i <= mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ for (i = 0, ap = mp->atoms, prp = mp->molprops; i <= mp->natoms + mp->nmolprops; i++) {
int s, t, ns, ne, mult;
Vector cr;
ne = s = t = 0;
if (i == mp->natoms) {
if (mp->cell == NULL || mp->frame_cells == NULL)
- break;
+ continue;
vp = mp->frame_cells;
mult = 4;
- } else {
+ } else if (i < mp->natoms) {
cr = ap->r;
vp = ap->frames;
mult = 1;
+ } else {
+ dp = prp->propvals;
}
for (j = 0; (ns = IntGroupGetStartPoint(group, j)) >= 0; j++) {
if (ns > ne) {
- memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (ns - ne));
+ if (i <= mp->natoms)
+ memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (ns - ne));
+ else
+ memmove((Double *)tempv + ne, dp + s, sizeof(Double) * (ns - ne));
s += ns - ne;
}
ne = IntGroupGetEndPoint(group, j);
tempv[ns * 4 + 2] = mp->cell->axes[2];
tempv[ns * 4 + 3] = mp->cell->origin;
}
- } else {
+ } else if (i < mp->natoms) {
if (inFrame != NULL)
tempv[ns] = inFrame[natoms * t + i];
else
tempv[ns] = cr;
+ } else {
+ ((Double *)tempv)[ns] = 0.0;
}
t++;
ns++;
}
}
if (n_new > ne) {
- memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (n_new - ne));
+ if (i <= mp->natoms)
+ memmove(tempv + ne * mult, vp + s * mult, sizeof(Vector) * mult * (n_new - ne));
+ else
+ memmove((Double *)tempv + ne, dp + s, sizeof(Double) * (n_new - ne));
s += n_new - ne;
}
if (i < mp->natoms)
ap->nframes = n_new;
- memmove(vp, tempv, sizeof(Vector) * mult * n_new);
+ if (i <= mp->natoms) {
+ memmove(vp, tempv, sizeof(Vector) * mult * n_new);
+ if (i < mp->natoms) {
+ ap->nframes = n_new;
+ ap = ATOM_NEXT(ap);
+ }
+ } else {
+ memmove(dp, (Double *)tempv, sizeof(Double) * n_new);
+ prp++;
+ }
}
free(tempv);
mp->nframes = n_new;
int i, count, n_new, n_old, natoms, nframes, old_count, new_cframe;
Vector *tempv, *vp;
Atom *ap;
+ MolProp *prp;
IntGroup *group, *group2;
if (mp == NULL || (natoms = mp->natoms) == 0 || (count = IntGroupGetCount(inGroup)) <= 0)
tempv[nl..n_old-1] -> ap[s..s+(n_old-nl-1)], s += n_old-nl
At last, s will become n_new and t will become count. */
nframes = 0;
- for (i = 0, ap = mp->atoms; i <= mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ for (i = 0, ap = mp->atoms, prp = mp->molprops; i <= mp->natoms + mp->nmolprops; i++) {
int s, t, j, ns, ne;
int mult;
/* if i == mp->natoms, mp->frame_cells is handled */
if (i == mp->natoms) {
if (mp->cell == NULL || mp->frame_cells == NULL)
- break;
- mult = 4;
+ continue;
+ mult = 4 * sizeof(Vector);
vp = mp->frame_cells;
old_count = n_old;
- } else {
- mult = 1;
+ } else if (i < mp->natoms) {
+ mult = sizeof(Vector);
vp = ap->frames;
if (vp == NULL) {
- ap->frames = vp = (Vector *)calloc(sizeof(Vector), n_old);
- if (vp == NULL) {
+ NewArray(&ap->frames, &ap->nframes, sizeof(Vector), n_old);
+ if (ap->frames == NULL) {
__MoleculeUnlock(mp);
return -1;
}
+ vp = ap->frames;
}
old_count = ap->nframes;
+ } else {
+ mult = sizeof(Double);
+ vp = (Vector *)prp->propvals;
+ old_count = n_old;
}
/* Copy vp to tempv */
- memset(tempv, 0, sizeof(Vector) * mult * n_old);
- memmove(tempv, vp, sizeof(Vector) * mult * (old_count > n_old ? n_old : old_count));
+ memset(tempv, 0, mult * n_old);
+ memmove(tempv, vp, mult * (old_count > n_old ? n_old : old_count));
ne = ns = s = t = 0;
for (j = 0; ns < n_old && (ns = IntGroupGetStartPoint(group, j)) >= 0; j++) {
if (ns > n_old)
ns = n_old;
if (ns > ne) {
- memmove(vp + s * mult, tempv + ne * mult, sizeof(Vector) * mult * (ns - ne));
+ memmove((char *)vp + s * mult, (char *)tempv + ne * mult, mult * (ns - ne));
s += ns - ne;
}
ne = IntGroupGetEndPoint(group, j);
while (ns < ne) {
if (i < mp->natoms)
outFrame[natoms * t + i] = tempv[ns];
- else if (outFrameCell != NULL) {
- outFrameCell[t * 4] = tempv[ns * 4];
- outFrameCell[t * 4 + 1] = tempv[ns * 4 + 1];
- outFrameCell[t * 4 + 2] = tempv[ns * 4 + 2];
- outFrameCell[t * 4 + 3] = tempv[ns * 4 + 3];
+ else if (i == mp->natoms) {
+ if (outFrameCell != NULL) {
+ outFrameCell[t * 4] = tempv[ns * 4];
+ outFrameCell[t * 4 + 1] = tempv[ns * 4 + 1];
+ outFrameCell[t * 4 + 2] = tempv[ns * 4 + 2];
+ outFrameCell[t * 4 + 3] = tempv[ns * 4 + 3];
+ }
}
t++;
ns++;
}
}
if (n_old > ne) {
- memmove(vp + s * mult, tempv + ne * mult, sizeof(Vector) * mult * (n_old - ne));
+ memmove((char *)vp + s * mult, (char *)tempv + ne * mult, mult * (n_old - ne));
s += n_old - ne;
}
if (i < mp->natoms)
free(ap->frames);
ap->frames = NULL;
ap->nframes = 0;
- } else {
+ } else if (i == mp->natoms) {
free(mp->frame_cells);
mp->frame_cells = NULL;
mp->nframe_cells = 0;
+ } else {
+ prp->propvals = (Double *)realloc(prp->propvals, sizeof(Double));
}
} else {
- if (i < mp->natoms)
- ap->frames = (Vector *)realloc(ap->frames, sizeof(Vector) * s);
- else {
+ if (i < mp->natoms) {
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), s - 1, NULL);
+ ap->nframes = s;
+ } else if (i == mp->natoms) {
AssignArray(&mp->frame_cells, &mp->nframe_cells, sizeof(Vector) * 4, s - 1, NULL);
mp->nframe_cells = s;
+ } else {
+ prp->propvals = (Double *)realloc(prp->propvals, sizeof(Double) * s);
}
}
+ if (i < mp->natoms) {
+ ap = ATOM_NEXT(ap);
+ } else if (i > mp->natoms) {
+ prp++;
+ }
}
free(tempv);
mp->nframes = nframes;
/* Write the current coordinate back to the frame array */
ap->frames[cframe] = ap->r;
}
- if (frame != cframe && frame >= 0 && frame < ap->nframes) {
+ if ((frame != cframe || copyback == 0) && frame >= 0 && frame < ap->nframes) {
/* Read the coordinate from the frame array */
ap->r = ap->frames[frame];
modified = 1;
vp[3] = mp->cell->origin;
}
/* Set the cell from the frame array */
- if (frame != cframe && frame >= 0 && frame < mp->nframe_cells) {
+ if ((frame != cframe || copyback == 0) && frame >= 0 && frame < mp->nframe_cells) {
MoleculeSetPeriodicBox(mp, &mp->frame_cells[frame * 4], &mp->frame_cells[frame * 4 + 1], &mp->frame_cells[frame * 4 + 2], &mp->frame_cells[frame * 4 + 3], mp->cell->flags, 0);
modified = 1;
MoleculeAmendBySymmetry(mp, NULL, NULL, NULL);
return nframes;
}
-#pragma mark ====== Pi Atoms ======
-
-#if PIATOM
int
-MoleculeCalculatePiAtomPosition(Molecule *mol, int idx)
+MoleculeReorderFrames(Molecule *mp, const Int *old_idx)
{
- PiAtom *pp;
- Int i, *cp;
- if (mol == NULL || idx < 0 || idx >= mol->npiatoms)
- return -1;
- pp = mol->piatoms + idx;
- cp = AtomConnectData(&pp->connect);
- VecZero(pp->r);
- for (i = pp->connect.count - 1; i >= 0; i--) {
- Vector rr = ATOM_AT_INDEX(mol->atoms, cp[i])->r;
- VecScaleInc(pp->r, rr, pp->coeffs[i]);
+ Int *ip, i, j, n, nframes;
+ Double *dp;
+ Atom *ap;
+ MolProp *prp;
+ if (mp == NULL || old_idx == NULL)
+ return 0;
+ nframes = MoleculeGetNumberOfFrames(mp);
+ MoleculeFlushFrames(mp);
+ ip = (Int *)malloc(sizeof(Int) * nframes);
+ if (ip == NULL)
+ return -1; /* Out of memory */
+ memset(ip, 0, sizeof(Int) * nframes);
+ /* Check the argument */
+ for (i = 0; i < nframes; i++) {
+ j = old_idx[i];
+ if (j < 0 || j >= nframes || ip[j] != 0) {
+ free(ip);
+ return -2; /* Bad argument */
+ }
+ ip[j] = 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;
- }
- }
+ free(ip);
+ dp = (Double *)malloc(sizeof(Double) * nframes * 12);
+ for (i = 0, ap = mp->atoms, prp = mp->molprops; i <= mp->natoms + mp->nmolprops; i++) {
+ for (j = 0; j < nframes; j++) {
+ n = old_idx[j];
+ if (i < mp->natoms) {
+ ((Vector *)dp)[j] = (n < ap->nframes ? ap->frames[n] : ap->r);
+ } else if (i == mp->natoms) {
+ if (mp->cell != NULL) {
+ if (n < mp->nframe_cells && mp->frame_cells != NULL)
+ memmove(dp + j * 12, mp->frame_cells + n * 4, sizeof(Vector) * 4);
+ else {
+ ((Vector *)dp)[j * 4] = mp->cell->axes[0];
+ ((Vector *)dp)[j * 4] = mp->cell->axes[1];
+ ((Vector *)dp)[j * 4] = mp->cell->axes[2];
+ ((Vector *)dp)[j * 4] = mp->cell->origin;
}
}
+ } else {
+ dp[j] = prp->propvals[n];
}
}
- 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;
+ for (j = 0; j < nframes; j++) {
+ if (i < mp->natoms) {
+ if (ap->nframes <= j)
+ AssignArray(&ap->frames, &ap->nframes, sizeof(Vector), nframes - 1, NULL);
+ ap->frames[j] = ((Vector *)dp)[j];
+ } else if (i == mp->natoms) {
+ if (mp->cell != NULL) {
+ AssignArray(&mp->frame_cells, &mp->nframe_cells, sizeof(Vector) * 4, nframes - 1, NULL);
+ memmove(mp->frame_cells + j * 4, dp + j * 12, sizeof(Vector) * 4);
+ }
+ } else {
+ prp->propvals[j] = dp[j];
}
- AssignArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), m - 1, NULL);
- for (j = mol->natoms + 1; j < m; j++)
- mol->piconnects[j] = -1;
}
+ if (i < mp->natoms)
+ ap = ATOM_NEXT(ap);
+ else if (i > mp->natoms)
+ prp++;
+ }
+ free(dp);
+ MoleculeSelectFrame(mp, mp->cframe, 0);
+ return 0;
+}
+
+#pragma mark ====== Molecule Propeties ======
+
+int
+MoleculeCreateProperty(Molecule *mp, const char *name)
+{
+ int i;
+ MolProp *prp;
+ for (i = 0, prp = mp->molprops; i < mp->nmolprops; i++, prp++) {
+ if (strcmp(prp->propname, name) == 0)
+ return -(i + 1);
+ }
+ prp = (MolProp *)calloc(sizeof(MolProp), 1);
+ if (prp == NULL)
+ return -10000;
+ prp->propname = strdup(name);
+ if (prp->propname == NULL)
+ return -10000;
+ i = MoleculeGetNumberOfFrames(mp);
+ prp->propvals = (Double *)calloc(sizeof(Double), i);
+ if (prp->propvals == NULL)
+ return -10000;
+ AssignArray(&mp->molprops, &mp->nmolprops, sizeof(MolProp), mp->nmolprops, prp);
+ free(prp);
+ return mp->nmolprops - 1;
+}
+
+int
+MoleculeLookUpProperty(Molecule *mp, const char *name)
+{
+ int i;
+ MolProp *prp;
+ for (i = 0, prp = mp->molprops; i < mp->nmolprops; i++, prp++) {
+ if (strcmp(prp->propname, name) == 0)
+ return i;
+ }
+ return -1;
+}
+
+int
+MoleculeDeletePropertyAtIndex(Molecule *mp, int idx)
+{
+ if (idx >= 0 && idx < mp->nmolprops) {
+ free(mp->molprops[idx].propname);
+ free(mp->molprops[idx].propvals);
+ DeleteArray(&mp->molprops, &mp->nmolprops, sizeof(MolProp), idx, 1, NULL);
+ return idx;
+ }
+ return -1;
+}
+
+int
+MoleculeSetProperty(Molecule *mp, int idx, IntGroup *ig, const Double *values)
+{
+ IntGroupIterator iter;
+ int i, n, nframes;
+ if (idx < 0 || idx >= mp->nmolprops)
+ return -1;
+ IntGroupIteratorInit(ig, &iter);
+ nframes = MoleculeGetNumberOfFrames(mp);
+ n = 0;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ if (i >= nframes)
+ break;
+ mp->molprops[idx].propvals[i] = values[n];
+ n++;
+ }
+ IntGroupIteratorRelease(&iter);
+ return n;
+}
+
+int
+MoleculeGetProperty(Molecule *mp, int idx, IntGroup *ig, Double *outValues)
+{
+ IntGroupIterator iter;
+ int i, n, nframes;
+ if (idx < 0 || idx >= mp->nmolprops)
+ return -1;
+ IntGroupIteratorInit(ig, &iter);
+ nframes = MoleculeGetNumberOfFrames(mp);
+ n = 0;
+ while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+ if (i >= nframes)
+ break;
+ outValues[n] = mp->molprops[idx].propvals[i];
+ n++;
}
- return (mol->npiconnects - mol->natoms - 1); /* Returns the number of entries */
+ IntGroupIteratorRelease(&iter);
+ return n;
+}
+
+#pragma mark ====== Pi Atoms ======
+
+static inline void
+sMoleculeCalculatePiAnchorPosition(Atom *ap, Atom *atoms)
+{
+ Int *cp, j, n;
+ Atom *ap2;
+ cp = AtomConnectData(&ap->anchor->connect);
+ n = ap->anchor->connect.count;
+ VecZero(ap->r);
+ for (j = 0; j < n; j++) {
+ Double w = ap->anchor->coeffs[j];
+ ap2 = ATOM_AT_INDEX(atoms, cp[j]);
+ VecScaleInc(ap->r, ap2->r, w);
+ }
}
void
-MoleculeInvalidatePiConnectionTable(Molecule *mol)
+MoleculeUpdatePiAnchorPositions(Molecule *mol)
{
- if (mol == NULL)
+ Int i;
+ Atom *ap;
+ for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
+ if (ap->anchor == NULL)
+ continue;
+ sMoleculeCalculatePiAnchorPosition(ap, mol->atoms);
+ }
+}
+
+void
+MoleculeCalculatePiAnchorPosition(Molecule *mol, int idx)
+{
+ Atom *ap;
+ if (mol == NULL || idx < 0 || idx >= mol->natoms)
+ return;
+ ap = ATOM_AT_INDEX(mol->atoms, idx);
+ if (ap->anchor == NULL)
return;
- if (mol->piconnects != NULL) {
- free(mol->piconnects);
- mol->piconnects = NULL;
+ sMoleculeCalculatePiAnchorPosition(ap, mol->atoms);
+}
+
+int
+MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Double *weights, Int *nUndoActions, struct MolAction ***undoActions)
+{
+ Atom *ap;
+ Int *ip, i, j, n, *np;
+ Double d;
+ if (mol == NULL || idx < 0 || idx >= mol->natoms || nentries <= 1)
+ return -1; /* Invalid argument */
+ if (weights != NULL) {
+ d = 0.0;
+ for (i = 0; i < nentries; i++) {
+ if (weights[i] <= 0.0) {
+ return 10; /* Weights must be positive */
+ }
+ d += weights[i];
+ }
+ d = 1.0 / d;
+ } else d = 1.0 / nentries;
+ ap = ATOM_AT_INDEX(mol->atoms, idx);
+ if (ap->anchor != NULL) {
+ /* Already an anchor: check if bonds/angles/dihedrals have entries related to this anchor */
+ IntGroup *bg, *ag, *dg, *ig;
+ Int *ibuf, ibufsize;
+ MolAction *act;
+ bg = ag = dg = ig = NULL;
+ ip = AtomConnectData(&ap->anchor->connect);
+ for (i = 0; i < ap->anchor->connect.count; i++) {
+ n = ip[i];
+ for (j = 0; j < nentries; j++) {
+ if (n == entries[j])
+ break;
+ }
+ if (j == nentries) {
+ /* This entry will disappear: if any bond/angle/dihedral has idx-n pair, that should be removed. */
+ for (j = 0, np = mol->bonds; j < mol->nbonds; j++, np += 2) {
+ if ((idx == np[0] && n == np[1]) || (idx == np[1] && n == np[0])) {
+ if (bg == NULL)
+ bg = IntGroupNew();
+ IntGroupAdd(bg, j, 1);
+ }
+ }
+ for (j = 0, np = mol->angles; j < mol->nangles; j++, np += 3) {
+ if ((idx == np[0] && n == np[1]) || (idx == np[1] && n == np[2]) ||
+ (idx == np[1] && n == np[0]) || (idx == np[2] && n == np[1])) {
+ if (ag == NULL)
+ ag = IntGroupNew();
+ IntGroupAdd(ag, j, 1);
+ }
+ }
+ for (j = 0, np = mol->dihedrals; j < mol->ndihedrals; j++, np += 4) {
+ if ((idx == np[0] && n == np[1]) || (idx == np[1] && n == np[2]) || (idx == np[2] && n == np[3]) ||
+ (idx == np[1] && n == np[0]) || (idx == np[2] && n == np[1]) || (idx == np[3] && n == np[2])) {
+ if (dg == NULL)
+ dg = IntGroupNew();
+ IntGroupAdd(dg, j, 1);
+ }
+ }
+ for (j = 0, np = mol->impropers; j < mol->nimpropers; j++, np += 4) {
+ if ((idx == np[0] && n == np[2]) || (idx == np[1] && n == np[2]) || (idx == np[3] && n == np[2]) ||
+ (idx == np[2] && n == np[0]) || (idx == np[2] && n == np[1]) || (idx == np[2] && n == np[3])) {
+ if (ig == NULL)
+ ig = IntGroupNew();
+ IntGroupAdd(ig, j, 1);
+ }
+ }
+ }
+ }
+ ibuf = NULL;
+ ibufsize = 0;
+ if (ig != NULL) {
+ /* Delete impropers (with undo info) */
+ i = IntGroupGetCount(ig);
+ AssignArray(&ibuf, &ibufsize, sizeof(Int), i * 4 - 1, NULL);
+ MoleculeDeleteImpropers(mol, ibuf, ig);
+ if (nUndoActions != NULL && undoActions != NULL) {
+ act = MolActionNew(gMolActionAddImpropers, i * 4, ibuf, ig);
+ AssignArray(undoActions, nUndoActions, sizeof(MolAction *), *nUndoActions, &act);
+ }
+ IntGroupRelease(ig);
+ }
+ if (dg != NULL) {
+ /* Delete dihedrals (with undo info) */
+ i = IntGroupGetCount(dg);
+ AssignArray(&ibuf, &ibufsize, sizeof(Int), i * 4 - 1, NULL);
+ MoleculeDeleteDihedrals(mol, ibuf, dg);
+ if (nUndoActions != NULL && undoActions != NULL) {
+ act = MolActionNew(gMolActionAddDihedrals, i * 4, ibuf, dg);
+ AssignArray(undoActions, nUndoActions, sizeof(MolAction *), *nUndoActions, &act);
+ }
+ IntGroupRelease(dg);
+ }
+ if (ag != NULL) {
+ /* Delete angles (with undo info) */
+ i = IntGroupGetCount(ag);
+ AssignArray(&ibuf, &ibufsize, sizeof(Int), i * 3 - 1, NULL);
+ MoleculeDeleteAngles(mol, ibuf, ag);
+ if (nUndoActions != NULL && undoActions != NULL) {
+ act = MolActionNew(gMolActionAddAngles, i * 3, ibuf, ag);
+ AssignArray(undoActions, nUndoActions, sizeof(MolAction *), *nUndoActions, &act);
+ }
+ IntGroupRelease(ag);
+ }
+ if (bg != NULL) {
+ /* Delete bonds (with undo info) */
+ i = IntGroupGetCount(bg);
+ AssignArray(&ibuf, &ibufsize, sizeof(Int), i * 2 - 1, NULL);
+ MoleculeDeleteBonds(mol, ibuf, bg, NULL, NULL);
+ if (nUndoActions != NULL && undoActions != NULL) {
+ act = MolActionNew(gMolActionAddBondsForUndo, i * 2, ibuf, bg);
+ AssignArray(undoActions, nUndoActions, sizeof(MolAction *), *nUndoActions, &act);
+ }
+ IntGroupRelease(bg);
+ }
+ } else {
+ ap->anchor = (PiAnchor *)calloc(sizeof(PiAnchor), 1);
+ }
+ AtomConnectResize(&ap->anchor->connect, nentries);
+ memmove(AtomConnectData(&ap->anchor->connect), entries, sizeof(Int) * nentries);
+ AssignArray(&ap->anchor->coeffs, &ap->anchor->ncoeffs, sizeof(Double), nentries - 1, NULL);
+ if (weights != NULL) {
+ memmove(ap->anchor->coeffs, weights, sizeof(Double) * nentries);
+ for (i = 0; i < nentries; i++)
+ ap->anchor->coeffs[i] *= d; /* Normalize weight */
+ } else {
+ for (i = 0; i < nentries; i++)
+ ap->anchor->coeffs[i] = d;
}
- mol->npiconnects = 0;
+ MoleculeCalculatePiAnchorPosition(mol, idx);
+ return 0;
}
-#endif /* PIATOM */
#pragma mark ====== MO calculation ======
/* Calculate an MO value for a single point. */
-/* Index is the MO number (1-based) */
+/* Index is the MO number (1-based); 0 denotes "arbitrary vector" */
/* tmp is an array of (natoms * 4) atoms, and used to store dr and |dr|^2 for each atom. */
static Double
-sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
+sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
{
ShellInfo *sp;
PrimInfo *pp;
Double val, tval, *cnp, *tmpp, *mobasep, *mop;
Int i, j;
/* Cache dr and |dr|^2 */
- for (i = 0; i < bset->natoms; i++) {
- Vector r = bset->pos[i];
- tmp[i * 4] = r.x = vp->x - r.x;
- tmp[i * 4 + 1] = r.y = vp->y - r.y;
- tmp[i * 4 + 2] = r.z = vp->z - r.z;
+ if (index == 0)
+ index = bset->nmos + 1;
+ for (i = 0; i < mp->natoms; i++) {
+ Vector r;
+ r = ATOM_AT_INDEX(mp->atoms, i)->r;
+ tmp[i * 4] = r.x = (vp->x - r.x) * kAngstrom2Bohr;
+ tmp[i * 4 + 1] = r.y = (vp->y - r.y) * kAngstrom2Bohr;
+ tmp[i * 4 + 2] = r.z = (vp->z - r.z) * kAngstrom2Bohr;
tmp[i * 4 + 3] = r.x * r.x + r.y * r.y + r.z * r.z;
}
/* Iterate over all shells */
for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
pp = bset->priminfos + sp->p_idx;
cnp = bset->cns + sp->cn_idx;
+ if (sp->a_idx >= mp->natoms)
+ return 0.0; /* This may happen when molecule is edited after setting up MO info */
tmpp = tmp + sp->a_idx * 4;
mop = mobasep + sp->m_idx;
switch (sp->sym) {
val += d0 + d1p + d1n + d2p + d2n;
break;
}
+ /* TODO: Support F/F7 and G/G9 type orbitals */
}
}
return val;
}
-/* Calculate one MO. The input vectors should be in bohr unit (angstrom * 1.889725989 = kAngstrom2Bohr). */
-/* mono is the MO number (1-based) */
+/* Calculate one MO. The input vectors are angstrom unit (changed from bohr unit: 20140520) */
+/* mono is the MO number (1-based); 0 denotes "arbitrary vector" */
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)
{
if (sSetupGaussianCoefficients(mp->bset) != 0)
return -1;
}
+ if (mp->bset->natoms_bs > mp->natoms)
+ return -3; /* Number of atoms is smaller than expected (internal error) */
+
cp = (Cube *)calloc(sizeof(Cube), 1);
if (cp == NULL) {
return -1;
cp->nz = nz;
/* TODO: use multithread */
- tmp = (Double *)calloc(sizeof(Double), mp->bset->natoms * 4);
+ tmp = (Double *)calloc(sizeof(Double), mp->bset->natoms_bs * 4);
if (tmp == NULL) {
free(cp->dp);
free(cp);
p.x = op->x + dxp->x * ix + dyp->x * iy + dzp->x * iz;
p.y = op->y + dxp->y * ix + dyp->y * iy + dzp->y * iz;
p.z = op->z + dxp->z * ix + dyp->z * iy + dzp->z * iz;
- cp->dp[n++] = sCalcMOPoint(mp->bset, mono, &p, tmp);
+ cp->dp[n++] = sCalcMOPoint(mp, mp->bset, mono, &p, tmp);
}
if (callback != NULL && n - nn > 100) {
nn = n;
return mp->bset->ncubes - 1;
}
+/* Output values are in angstrom unit (changed from bohr unit: 20140520) */
int
MoleculeGetDefaultMOGrid(Molecule *mp, Int npoints, Vector *op, Vector *xp, Vector *yp, Vector *zp, Int *nx, Int *ny, Int *nz)
{
int i;
- Vector rmin, rmax, *vp;
+ Vector rmin, rmax, r;
Double dr, dx, dy, dz;
- if (mp == NULL || mp->bset == NULL || mp->bset->natoms == 0)
+ Atom *ap;
+ if (mp == NULL || mp->bset == NULL)
return -1;
if (npoints <= 0)
npoints = 1000000;
rmin.x = rmin.y = rmin.z = 1e10;
rmax.x = rmax.y = rmax.z = -1e10;
- for (i = 0, vp = mp->bset->pos; i < mp->bset->natoms; i++, vp++) {
- dr = RadiusForAtomicNumber(ATOM_AT_INDEX(mp->atoms, i)->atomicNumber);
+ for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+ dr = RadiusForAtomicNumber(ap->atomicNumber);
+ r = ap->r;
if (dr == 0.0)
dr = 1.0;
- dr = dr * kAngstrom2Bohr * 3.0 + 2.0;
- if (rmin.x > vp->x - dr)
- rmin.x = vp->x - dr;
- if (rmin.y > vp->y - dr)
- rmin.y = vp->y - dr;
- if (rmin.z > vp->z - dr)
- rmin.z = vp->z - dr;
- if (rmax.x < vp->x + dr)
- rmax.x = vp->x + dr;
- if (rmax.y < vp->y + dr)
- rmax.y = vp->y + dr;
- if (rmax.z < vp->z + dr)
- rmax.z = vp->z + dr;
+ dr = dr * 3.0 + 2.0;
+ if (rmin.x > r.x - dr)
+ rmin.x = r.x - dr;
+ if (rmin.y > r.y - dr)
+ rmin.y = r.y - dr;
+ if (rmin.z > r.z - dr)
+ rmin.z = r.z - dr;
+ if (rmax.x < r.x + dr)
+ rmax.x = r.x + dr;
+ if (rmax.y < r.y + dr)
+ rmax.y = r.y + dr;
+ if (rmax.z < r.z + dr)
+ rmax.z = r.z + dr;
}
dx = rmax.x - rmin.x;
dy = rmax.y - rmin.y;
fprintf(fp, "%s MO=%d\n", comment, cp->idn);
fprintf(fp, " MO coefficients\n");
- fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", -(mp->bset->natoms), cp->origin.x, cp->origin.y, cp->origin.z);
- fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->nx, cp->dx.x, cp->dx.y, cp->dx.z);
- fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->ny, cp->dy.x, cp->dy.y, cp->dy.z);
- fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->nz, cp->dz.x, cp->dz.y, cp->dz.z);
+ fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", -(mp->bset->natoms_bs),
+ cp->origin.x * kAngstrom2Bohr, cp->origin.y * kAngstrom2Bohr, cp->origin.z * kAngstrom2Bohr);
+ fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->nx,
+ cp->dx.x * kAngstrom2Bohr, cp->dx.y * kAngstrom2Bohr, cp->dx.z * kAngstrom2Bohr);
+ fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->ny,
+ cp->dy.x * kAngstrom2Bohr, cp->dy.y * kAngstrom2Bohr, cp->dy.z * kAngstrom2Bohr);
+ fprintf(fp, "%5d %11.6f %11.6f %11.6f\n", cp->nz,
+ cp->dz.x * kAngstrom2Bohr, cp->dz.y * kAngstrom2Bohr, cp->dz.z * kAngstrom2Bohr);
/* Atomic information */
- for (i = 0; i < mp->bset->natoms; i++) {
- Vector *vp = mp->bset->pos + i;
+ for (i = 0; i < mp->natoms; i++) {
Atom *ap = ATOM_AT_INDEX(mp->atoms, i);
/* The second number should actually be the effective charge */
- fprintf(fp, "%5d %11.6f %11.6f %11.6f %11.6f\n", ap->atomicNumber, (double)ap->atomicNumber, vp->x, vp->y, vp->z);
+ fprintf(fp, "%5d %11.6f %11.6f %11.6f %11.6f\n", ap->atomicNumber, (double)ap->atomicNumber,
+ ap->r.x * kAngstrom2Bohr, ap->r.y * kAngstrom2Bohr, ap->r.z * kAngstrom2Bohr);
}
fprintf(fp, "%5d%5d\n", 1, 1);
for (i = n = 0; i < cp->nx; i++) {
for (j = 0; j < cp->ny; j++) {
for (k = 0; k < cp->nz; k++) {
- fprintf(fp, " %12.5e", cp->dp[n++]);
+ /* On Windows, the "%e" format writes the exponent in 3 digits, but
+ this is not standard. So we avoid using %e */
+ Double d = cp->dp[n++];
+ int exponent;
+ Double base;
+ if (d >= -1.0e-90 && d <= 1.0e-90) {
+ exponent = 0;
+ base = 0.0;
+ } else {
+ exponent = (int)floor(log10(fabs(d)));
+ base = d * pow(10, -1.0 * exponent);
+ }
+ fprintf(fp, " %8.5fe%+03d", base, exponent);
+ /* fprintf(fp, " %12.5e", d); */
if (k == cp->nz - 1 || k % 6 == 5)
fprintf(fp, "\n");
}
fclose(fp);
return 0;
}
+
+#pragma mark ====== Marching Cube (for isosurface) ======
+
+MCube *
+MoleculeClearMCube(Molecule *mol, Int nx, Int ny, Int nz, const Vector *origin, Double dx, Double dy, Double dz)
+{
+ MCube *mc = mol->mcube;
+ int i;
+ float rgba[8] = { 1, 1, 1, 0.6, 0, 0, 1, 0.6 };
+ if (mc != NULL) {
+ free(mc->dp);
+ free(mc->radii);
+ free(mc->c[0].fp);
+ free(mc->c[0].cubepoints);
+ free(mc->c[0].triangles);
+ free(mc->c[1].fp);
+ free(mc->c[1].cubepoints);
+ free(mc->c[1].triangles);
+ memmove(rgba, mc->c[0].rgba, sizeof(float) * 4);
+ memmove(rgba + 4, mc->c[1].rgba, sizeof(float) * 4);
+ free(mc);
+ mol->mcube = NULL;
+ }
+ if (nx > 0 && ny > 0 && nz > 0) {
+ mc = (MCube *)calloc(sizeof(MCube), 1);
+ mc->idn = -1;
+ /* round up to nearest 4N+1 integer */
+ dx *= nx;
+ dy *= ny;
+ dz *= nz;
+ mc->nx = (nx + 2) / 4 * 4 + 1;
+ mc->ny = (ny + 2) / 4 * 4 + 1;
+ mc->nz = (nz + 2) / 4 * 4 + 1;
+ mc->dx = dx / mc->nx;
+ mc->dy = dy / mc->ny;
+ mc->dz = dz / mc->nz;
+ mc->origin = *origin;
+ mc->dp = (Double *)malloc(sizeof(Double) * mc->nx * mc->ny * mc->nz);
+ if (mc->dp == NULL) {
+ free(mc);
+ return NULL;
+ }
+ mc->radii = (Double *)calloc(sizeof(Double), mol->natoms);
+ if (mc->radii == NULL) {
+ free(mc->dp);
+ free(mc);
+ return NULL;
+ }
+ mc->nradii = mol->natoms;
+ mc->c[0].fp = (unsigned char *)calloc(sizeof(unsigned char), mc->nx * mc->ny * mc->nz);
+ mc->c[1].fp = (unsigned char *)calloc(sizeof(unsigned char), mc->nx * mc->ny * mc->nz);
+ if (mc->c[0].fp == NULL || mc->c[1].fp == NULL) {
+ free(mc->c[0].fp);
+ free(mc->c[1].fp);
+ free(mc->dp);
+ free(mc->radii);
+ free(mc);
+ return NULL;
+ }
+ for (i = 0; i < mc->nx * mc->ny * mc->nz; i++) {
+ mc->dp[i] = DBL_MAX;
+ }
+ memmove(mc->c[0].rgba, rgba, sizeof(float) * 4);
+ memmove(mc->c[1].rgba, rgba + 4, sizeof(float) * 4);
+ mol->mcube = mc;
+ }
+ MoleculeCallback_notifyModification(mol, 0);
+ return mol->mcube;
+}
+
+static int sMarchingCubeTable[256][16] = {
+ {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+ {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+ {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+ {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
+ {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+ {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+ {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+ {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
+ {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
+ {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
+ {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+ {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+ {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+ {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
+ {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+ {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+ {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
+ {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
+ {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+ {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+ {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+ {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
+ {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
+ {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
+ {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+ {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+ {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
+ {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
+ {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
+ {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
+ {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
+ {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
+ {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
+ {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+ {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
+ {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+ {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+ {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
+ {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
+ {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
+ {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+ {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+ {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
+ {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
+ {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+ {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
+ {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
+ {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
+ {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
+ {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+ {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+ {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
+ {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+ {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
+ {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+ {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
+ {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
+ {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
+ {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+ {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+ {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
+ {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
+ {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+ {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
+ {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+ {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+ {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
+ {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
+ {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
+ {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
+ {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
+ {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+ {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+ {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+ {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
+ {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+ {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
+ {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
+ {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+ {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
+ {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
+ {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+ {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+ {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
+ {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
+ {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
+ {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
+ {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+ {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
+ {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+ {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
+ {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+ {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
+ {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+ {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+ {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
+ {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+ {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
+ {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
+ {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
+ {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
+ {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
+ {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
+ {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
+ {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
+ {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
+ {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
+ {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+ {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
+ {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
+ {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
+ {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
+ {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
+ {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
+ {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
+ {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+ {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
+ {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
+ {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
+ {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
+ {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
+ {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
+ {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
+ {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
+ {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
+ {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+ {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
+ {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
+ {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
+ {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+ {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+ {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+ {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
+ {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
+ {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
+ {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
+ {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
+ {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
+ {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
+ {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+ {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
+ {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
+ {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
+ {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+ {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
+ {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
+ {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
+ {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
+ {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
+ {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+ {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
+ {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
+ {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
+ {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
+ {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+ {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
+ {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+ {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+ {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+ {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
+ {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+ {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
+ {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+ {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
+};
+
+/* Recalculate the MCube */
+/* If idn < 0, then the current grid settings and values are unchanged, and */
+/* only the marching cubes are regenerated. */
+int
+MoleculeUpdateMCube(Molecule *mol, int idn)
+{
+ Int retval, step, sn;
+ Int n, ix, iy, iz, nx, ny, nz;
+ Int nn, iix, iiy, iiz;
+ Int ncubepoints, c1, c2, c3;
+ Int *ip;
+ Double thres, *tmp, dd;
+ Vector p;
+ MCube *mc;
+ MCubePoint *mcp;
+ Atom *ap;
+
+ if (mol == NULL || mol->bset == NULL || mol->mcube == NULL)
+ return -1;
+ if (mol->bset->cns == NULL) {
+ if (sSetupGaussianCoefficients(mol->bset) != 0)
+ return -1;
+ }
+ if (mol->bset->natoms_bs > mol->natoms)
+ return -1; /* Number of atoms is smaller than expected */
+
+ mc = mol->mcube;
+ if (idn >= 0) {
+ ShellInfo *sp;
+ Double *mobasep, *mop, mopmax;
+ Double xmin, xmax, ymin, ymax, zmin, zmax;
+ /* Clear mcube values */
+ for (ix = 0; ix < mc->nx * mc->ny * mc->nz; ix++) {
+ mc->dp[ix] = DBL_MAX;
+ mc->c[0].fp[ix] = 0;
+ mc->c[1].fp[ix] = 0;
+ }
+ mc->idn = idn;
+ /* Estimate the orbital sizes */
+ mc->radii = (Double *)realloc(mc->radii, sizeof(Double) * mol->natoms);
+ if (mc->radii == NULL)
+ return -2; /* Out of memory */
+ mc->nradii = mol->natoms;
+ if (mc->idn == mol->bset->nmos + 1) {
+ /* Total electron density */
+ for (ix = 0; ix < mol->natoms; ix++)
+ mc->radii[ix] = 1.0;
+ mopmax = 1.0;
+ } else {
+ memset(mc->radii, 0, sizeof(Double) * mc->nradii);
+ mobasep = mol->bset->mo + (mc->idn == 0 ? mol->bset->nmos : mc->idn - 1) * mol->bset->ncomps;
+ mopmax = 0.0;
+ for (ix = 0, sp = mol->bset->shells; ix < mol->bset->nshells; ix++, sp++) {
+ if (sp->a_idx >= mol->natoms)
+ continue; /* This may happen when molecule is edited after setting up MO info */
+ mop = mobasep + sp->m_idx;
+ for (iy = 0; iy < sp->ncomp; iy++) {
+ dd = fabs(mop[iy]);
+ if (dd > mc->radii[sp->a_idx])
+ mc->radii[sp->a_idx] = dd;
+ if (dd > mopmax)
+ mopmax = dd;
+ }
+ }
+ }
+ xmin = ymin = zmin = 1e10;
+ xmax = ymax = zmax = -1e10;
+ for (ix = 0, ap = mol->atoms; ix < mol->natoms; ix++, ap = ATOM_NEXT(ap)) {
+ dd = RadiusForAtomicNumber(ap->atomicNumber);
+ dd = (dd * 2.0 + 1.0) * (mc->radii[ix] / mopmax) * (mc->expand > 0.0 ? mc->expand : 1.0);
+ mc->radii[ix] = dd;
+ p = ap->r;
+ dd += 0.1;
+ if (p.x - dd < xmin)
+ xmin = p.x - dd;
+ if (p.y - dd < ymin)
+ ymin = p.y - dd;
+ if (p.z - dd < zmin)
+ zmin = p.z - dd;
+ if (p.x + dd > xmax)
+ xmax = p.x + dd;
+ if (p.y + dd > ymax)
+ ymax = p.y + dd;
+ if (p.z + dd > zmax)
+ zmax = p.z + dd;
+ }
+ mc->origin.x = xmin;
+ mc->origin.y = ymin;
+ mc->origin.z = zmin;
+ mc->dx = (xmax - xmin) / mc->nx;
+ mc->dy = (ymax - ymin) / mc->ny;
+ mc->dz = (zmax - zmin) / mc->nz;
+ }
+
+ /* Temporary work area */
+ tmp = (Double *)calloc(sizeof(Double), mol->bset->natoms_bs * 4);
+ if (tmp == NULL)
+ return -2;
+
+ /* TODO: use multithread */
+ nx = mc->nx;
+ ny = mc->ny;
+ nz = mc->nz;
+ step = 4;
+
+#if 1
+ /* Calculate points within certain distances from atoms */
+ for (nn = 0, ap = mol->atoms; nn < mol->natoms; nn++, ap = ATOM_NEXT(ap)) {
+ /* dd = RadiusForAtomicNumber(ap->atomicNumber);
+ if (dd == 0.0)
+ dd = 1.0;
+ dd = dd * 1.5 + 1.0; */
+ dd = mc->radii[nn];
+ p.x = ap->r.x - dd - mc->origin.x;
+ p.y = ap->r.y - dd - mc->origin.y;
+ p.z = ap->r.z - dd - mc->origin.z;
+ c1 = p.x / mc->dx;
+ c2 = p.y / mc->dy;
+ c3 = p.z / mc->dz;
+ iix = c1 + ceil(dd * 2.0 / mc->dx);
+ iiy = c2 + ceil(dd * 2.0 / mc->dy);
+ iiz = c3 + ceil(dd * 2.0 / mc->dz);
+ if (c1 < 0)
+ c1 = 0;
+ if (c2 < 0)
+ c2 = 0;
+ if (c3 < 0)
+ c3 = 0;
+ if (iix >= nx)
+ iix = nx - 1;
+ if (iiy >= ny)
+ iiy = ny - 1;
+ if (iiz >= nz)
+ iiz = nz - 1;
+ for (ix = c1; ix <= iix; ix++) {
+ p.x = mc->origin.x + mc->dx * ix;
+ for (iy = c2; iy <= iiy; iy++) {
+ p.y = mc->origin.y + mc->dy * iy;
+ for (iz = c3; iz <= iiz; iz++) {
+ n = (ix * ny + iy) * nz + iz;
+ if (mc->dp[n] == DBL_MAX) {
+ p.z = mc->origin.z + mc->dz * iz;
+ if (mc->idn == mol->bset->nmos + 1) {
+ /* Total electron density */
+ Int ne_alpha, ne_beta;
+ mc->dp[n] = 0.0;
+ ne_alpha = mol->bset->ne_alpha;
+ ne_beta = mol->bset->ne_beta;
+ if (mol->bset->rflag == 2 && ne_alpha < ne_beta) {
+ /* ROHF case: ensure ne_alpha >= ne_beta */
+ ne_beta = ne_alpha;
+ ne_alpha = mol->bset->ne_beta;
+ }
+ for (sn = 1; sn <= ne_alpha; sn++) {
+ dd = sCalcMOPoint(mol, mol->bset, sn, &p, tmp);
+ dd = dd * dd;
+ if (mol->bset->rflag != 0 && sn <= ne_beta)
+ dd *= 2;
+ mc->dp[n] += dd;
+ }
+ if (mol->bset->rflag == 0) {
+ for (sn = 1; sn <= ne_beta; sn++) {
+ dd = sCalcMOPoint(mol, mol->bset, sn + mol->bset->ncomps, &p, tmp);
+ mc->dp[n] += dd * dd;
+ }
+ }
+ } else {
+ mc->dp[n] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+ }
+ }
+ }
+ }
+ }
+ }
+
+#else
+ /* (i * step, j * step, k * step) */
+ for (ix = 0; ix < nx; ix += step) {
+ for (iy = 0; iy < ny; iy += step) {
+ for (iz = 0; iz < nz; iz += step) {
+ n = (ix * ny + iy) * nz + iz;
+ if (mc->dp[n] == DBL_MAX) {
+ p.x = mc->origin.x + mc->dx * ix;
+ p.y = mc->origin.y + mc->dy * iy;
+ p.z = mc->origin.z + mc->dz * iz;
+ mc->dp[n] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+ }
+ n += step;
+ }
+ }
+ }
+
+ /* Intermediate points */
+ for (step = 4; step > 1; step /= 2) {
+ hstep = step / 2;
+ for (sn = 0; sn <= 1; sn++) {
+ n = 0;
+ for (ix = 0; ix < nx - 1; ix += step) {
+ for (iy = 0; iy < ny - 1; iy += step) {
+ for (iz = 0; iz < nz - 1; iz += step) {
+ flags = 0;
+ thres = mc->thres * (sn == 0 ? 1 : -1);
+ n = (ix * ny + iy) * nz + iz;
+ if (mc->dp[n] == DBL_MAX || mc->dp[n + step * (nz * (ny + 1) + 1)] == DBL_MAX)
+ continue;
+ /* (ix, iy, iz) */
+ if (mc->dp[n] >= thres)
+ flags |= 1;
+ /* (ix + step, iy, iz) */
+ if (mc->dp[n + step * ny * nz] >= thres)
+ flags |= 2;
+ /* (ix, iy + step, iz) */
+ if (mc->dp[n + step * nz] >= thres)
+ flags |= 4;
+ /* (ix + 4, iy + step, iz) */
+ if (mc->dp[n + step * nz * (ny + 1)] >= thres)
+ flags |= 8;
+ /* (ix, iy, iz + step) */
+ if (mc->dp[n + step] >= thres)
+ flags |= 16;
+ if (mc->dp[n + step * (ny * nz + 1)] >= thres)
+ flags |= 32;
+ /* (ix, iy + step, iz + step) */
+ if (mc->dp[n + step * (nz + 1)] >= thres)
+ flags |= 64;
+ /* (ix + step, iy + step, iz + step) */
+ if (mc->dp[n + step * (nz * (ny + 1) + 1)] >= thres)
+ flags |= 128;
+ if (flags != 0 && flags != 255) {
+ /* Calc the intermediate points */
+ for (iix = 0; iix <= step; iix += hstep) {
+ for (iiy = 0; iiy <= step; iiy += hstep) {
+ for (iiz = 0; iiz <= step; iiz += hstep) {
+ if (iix % step == 0 && iiy % step == 0 && iiz % step == 0)
+ continue;
+ nn = n + (iix * ny + iiy) * nz + iiz;
+ if (mc->dp[nn] == DBL_MAX) {
+ p.x = mc->origin.x + mc->dx * (ix + iix);
+ p.y = mc->origin.y + mc->dy * (iy + iiy);
+ p.z = mc->origin.z + mc->dz * (iz + iiz);
+ mc->dp[nn] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+#endif
+
+ free(tmp);
+
+ /* Calculate vertex positions and normal vectors */
+ for (sn = 0; sn <= 1; sn++) {
+ n = 0;
+ thres = mc->thres * (sn == 0 ? 1 : -1);
+ VecZero(p);
+ for (ix = 0; ix < nx - 1; ix++) {
+ for (iy = 0; iy < ny - 1; iy++) {
+ for (iz = 0; iz < nz - 1; iz++) {
+ Double dd0, dd1;
+ nn = (ix * ny + iy) * nz + iz;
+ dd0 = mc->dp[nn];
+ if (dd0 == DBL_MAX)
+ continue;
+ if (0) {
+ dd1 = mc->dp[nn + ny * nz];
+ if (dd1 != DBL_MAX)
+ p.x = (dd1 - dd0) / mc->dx;
+ else if (ix > 0 && (dd1 = mc->dp[nn - ny * nz]) != DBL_MAX)
+ p.x = (dd0 - dd1) / mc->dx;
+ else continue; /* Cannot define gradient */
+ dd1 = mc->dp[nn + nz];
+ if (dd1 != DBL_MAX)
+ p.y = (dd1 - dd0) / mc->dy;
+ else if (iy > 0 && (dd1 = mc->dp[nn - nz]) != DBL_MAX)
+ p.y = (dd0 - dd1) / mc->dy;
+ else continue;
+ dd1 = mc->dp[nn + 1];
+ if (dd1 != DBL_MAX)
+ p.z = (dd1 - dd0) / mc->dz;
+ else if (iz > 0 && (dd1 = mc->dp[nn - 1]) != DBL_MAX)
+ p.z = (dd0 - dd1) / mc->dz;
+ else continue;
+ NormalizeVec(&p, &p);
+ }
+ if (n + 3 >= mc->c[sn].ncubepoints) {
+ /* Expand cubepoints[] array */
+ mc->c[sn].cubepoints = (MCubePoint *)realloc(mc->c[sn].cubepoints, sizeof(MCubePoint) * (mc->c[sn].ncubepoints + 8192));
+ if (mc->c[sn].cubepoints == NULL) {
+ mc->c[sn].ncubepoints = 0;
+ retval = -3;
+ goto end;
+ }
+ mc->c[sn].ncubepoints += 8192;
+ }
+ mcp = mc->c[sn].cubepoints + n;
+ iix = (dd0 >= thres ? 1 : -1);
+ /* (x, y, z)->(x + 1, y, z) */
+ dd1 = mc->dp[nn + ny * nz];
+ if (dd1 != DBL_MAX) {
+ iiy = (dd1 >= thres ? 1 : -1);
+ if (iix != iiy) {
+ /* Register */
+ mcp->key = nn * 3;
+ mcp->d = (thres - dd0) / (dd1 - dd0);
+ mcp->pos[0] = mc->origin.x + mc->dx * (ix + mcp->d);
+ mcp->pos[1] = mc->origin.y + mc->dy * iy;
+ mcp->pos[2] = mc->origin.z + mc->dz * iz;
+ mcp->grad[0] = p.x;
+ mcp->grad[1] = p.y;
+ mcp->grad[2] = p.z;
+ mcp++;
+ n++;
+ }
+ }
+ /* (x, y, z)->(x, y + 1, z) */
+ dd1 = mc->dp[nn + nz];
+ if (dd1 != DBL_MAX) {
+ iiy = (dd1 >= thres ? 1 : -1);
+ if (iix != iiy) {
+ /* Register */
+ mcp->key = nn * 3 + 1;
+ mcp->d = (thres - dd0) / (dd1 - dd0);
+ mcp->pos[0] = mc->origin.x + mc->dx * ix;
+ mcp->pos[1] = mc->origin.y + mc->dy * (iy + mcp->d);
+ mcp->pos[2] = mc->origin.z + mc->dz * iz;
+ mcp->grad[0] = p.x;
+ mcp->grad[1] = p.y;
+ mcp->grad[2] = p.z;
+ mcp++;
+ n++;
+ }
+ }
+ /* (x, y, z)->(x, y, z + 1) */
+ dd1 = mc->dp[nn + 1];
+ if (dd1 != DBL_MAX) {
+ iiy = (dd1 >= thres ? 1 : -1);
+ if (iix != iiy) {
+ /* Register */
+ mcp->key = nn * 3 + 2;
+ mcp->d = (thres - dd0) / (dd1 - dd0);
+ mcp->pos[0] = mc->origin.x + mc->dx * ix;
+ mcp->pos[1] = mc->origin.y + mc->dy * iy;
+ mcp->pos[2] = mc->origin.z + mc->dz * (iz + mcp->d);
+ mcp->grad[0] = p.x;
+ mcp->grad[1] = p.y;
+ mcp->grad[2] = p.z;
+ mcp++;
+ n++;
+ }
+ }
+ }
+ }
+ }
+ if (n < mc->c[sn].ncubepoints)
+ mc->c[sn].cubepoints[n].key = -1; /* End mark */
+ ncubepoints = n;
+ if (ncubepoints < 3) {
+ /* Less than 3 points: no triangles */
+ if (mc->c[sn].ntriangles > 0)
+ mc->c[sn].triangles[0] = -1; /* End mark */
+ continue;
+ }
+
+ /* Create triangle table */
+ n = 0;
+ for (ix = 0; ix < nx - 1; ix++) {
+ for (iy = 0; iy < ny - 1; iy++) {
+ for (iz = 0; iz < nz - 1; iz++) {
+ nn = (ix * ny + iy) * nz + iz;
+ iix = 0;
+ if ((dd = mc->dp[nn]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 1;
+ if ((dd = mc->dp[nn + ny * nz]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 2;
+ if ((dd = mc->dp[nn + ny * nz + nz]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 4;
+ if ((dd = mc->dp[nn + nz]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 8;
+ if ((dd = mc->dp[nn + 1]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 16;
+ if ((dd = mc->dp[nn + ny * nz + 1]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 32;
+ if ((dd = mc->dp[nn + ny * nz + nz + 1]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 64;
+ if ((dd = mc->dp[nn + nz + 1]) == DBL_MAX)
+ continue;
+ else if (dd >= thres)
+ iix |= 128;
+ for (iiy = 0; iiy < 15; iiy++) {
+ nn = sMarchingCubeTable[iix][iiy];
+ if (nn < 0)
+ break;
+ /* key index for edges 0-11 */
+ switch (nn) {
+ case 0: iiz = (( ix * ny + iy ) * nz + iz ) * 3; break;
+ case 1: iiz = (((ix + 1) * ny + iy ) * nz + iz ) * 3 + 1; break;
+ case 2: iiz = (( ix * ny + iy + 1) * nz + iz ) * 3; break;
+ case 3: iiz = (( ix * ny + iy ) * nz + iz ) * 3 + 1; break;
+ case 4: iiz = (( ix * ny + iy ) * nz + iz + 1) * 3; break;
+ case 5: iiz = (((ix + 1) * ny + iy ) * nz + iz + 1) * 3 + 1; break;
+ case 6: iiz = (( ix * ny + iy + 1) * nz + iz + 1) * 3; break;
+ case 7: iiz = (( ix * ny + iy ) * nz + iz + 1) * 3 + 1; break;
+ case 8: iiz = (( ix * ny + iy ) * nz + iz ) * 3 + 2; break;
+ case 9: iiz = (((ix + 1) * ny + iy ) * nz + iz ) * 3 + 2; break;
+ case 10: iiz = (((ix + 1) * ny + iy + 1) * nz + iz ) * 3 + 2; break;
+ case 11: iiz = (( ix * ny + iy + 1) * nz + iz ) * 3 + 2; break;
+ default:
+ /* Skip this triangle */
+ iiy = (iiy - iiy % 3) + 2;
+ n = n - n % 3;
+ continue;
+ }
+ /* Look for the key index in cubepoints */
+ c1 = 0;
+ c3 = ncubepoints - 1;
+ mcp = mc->c[sn].cubepoints;
+ while (1) {
+ int w;
+ /* c1 is always less than c3 */
+ if (c1 + 1 == c3) {
+ /* end of search */
+ if (mcp[c1].key == iiz) {
+ c2 = c1;
+ } else if (mcp[c3].key == iiz) {
+ c2 = c3;
+ } else {
+ c2 = -1;
+ }
+ break;
+ }
+ c2 = (c1 + c3) / 2;
+ w = mcp[c2].key - iiz;
+ if (w == 0)
+ break;
+ if (w < 0) {
+ c1 = c2;
+ } else {
+ c3 = c2;
+ }
+ }
+ if (c2 < 0) {
+ /* Not found: skip this triangle */
+ iiy = (iiy - iiy % 3) + 2;
+ n = n - n % 3;
+ continue;
+ }
+ if (n + 1 >= mc->c[sn].ntriangles) {
+ /* Expand triangles[] array */
+ mc->c[sn].triangles = (Int *)realloc(mc->c[sn].triangles, sizeof(Int) * (mc->c[sn].ntriangles + 8192));
+ if (mc->c[sn].triangles == NULL) {
+ mc->c[sn].ntriangles = 0;
+ retval = -4;
+ goto end;
+ }
+ mc->c[sn].ntriangles += 8192;
+ }
+ mc->c[sn].triangles[n] = c2;
+ n++;
+ }
+ }
+ }
+ }
+ if (n < mc->c[sn].ntriangles)
+ mc->c[sn].triangles[n] = -1; /* End mark */
+
+ /* Estimate the normal vector */
+ for (n = 0, ip = mc->c[sn].triangles; ip[n] >= 0; n += 3) {
+ Vector v[3];
+ for (ix = 0; ix < 3; ix++) {
+ mcp = &(mc->c[sn].cubepoints[ip[n + ix]]);
+ v[ix].x = mcp->pos[0];
+ v[ix].y = mcp->pos[1];
+ v[ix].z = mcp->pos[2];
+ }
+ VecDec(v[2], v[0]);
+ VecDec(v[1], v[0]);
+ VecCross(v[0], v[1], v[2]);
+ NormalizeVec(v, v);
+ for (ix = 0; ix < 3; ix++) {
+ mcp = &(mc->c[sn].cubepoints[ip[n + ix]]);
+ mcp->grad[0] += v[0].x;
+ mcp->grad[1] += v[0].y;
+ mcp->grad[2] += v[0].z;
+ }
+ }
+ for (n = 0, mcp = mc->c[sn].cubepoints; mcp->key >= 0; mcp++) {
+ if (mcp->grad[0] != 0.0 || mcp->grad[1] != 0.0 || mcp->grad[2] != 0.0) {
+ dd = 1.0 / sqrt(mcp->grad[0] * mcp->grad[0] + mcp->grad[1] * mcp->grad[1] + mcp->grad[2] * mcp->grad[2]);
+ if (mc->thres < 0.0)
+ dd = -dd;
+ mcp->grad[0] *= dd;
+ mcp->grad[1] *= dd;
+ mcp->grad[2] *= dd;
+ }
+ }
+ }
+ retval = 0;
+ MoleculeCallback_notifyModification(mol, 0);
+end:
+ /* For debug */
+ if (0) {
+ char *MyAppCallback_getDocumentHomeDir(void);
+ FILE *fp;
+ char *s;
+ Double dmax, dmin;
+ asprintf(&s, "%s/%s", MyAppCallback_getDocumentHomeDir(), "mcube_log.txt");
+ fp = fopen(s, "w");
+ dmax = -1e8;
+ dmin = 1e8;
+ for (n = 0; n < mc->nx * mc->ny * mc->nz; n++) {
+ if (mc->dp[n] == DBL_MAX)
+ continue;
+ if (dmax < mc->dp[n])
+ dmax = mc->dp[n];
+ if (dmin > mc->dp[n])
+ dmin = mc->dp[n];
+ }
+ dmax = fabs(dmax);
+ dmin = fabs(dmin);
+ if (dmax < dmin)
+ dmax = dmin;
+ dmax = 1.001 * dmax;
+ fprintf(fp, "thres = %g = 100\n", mc->thres);
+ for (iz = 0; iz < mc->nz; iz++) {
+ fprintf(fp, "z = %d\n", iz);
+ for (iy = 0; iy < mc->ny; iy++) {
+ for (ix = 0; ix < mc->nx; ix++) {
+ n = (ix * ny + iy) * nz + iz;
+ dd = mc->dp[n];
+ if (dd == DBL_MAX)
+ fprintf(fp, " XXX ");
+ else {
+ dd = dd * 100 / mc->thres;
+ if (dd > 999.0)
+ dd = 999.0;
+ else if (dd < -999.0)
+ dd = -999.0;
+ fprintf(fp, "%4d ", (int)(dd));
+ }
+ }
+ fprintf(fp, "\n");
+ }
+ fprintf(fp, "\n");
+ }
+
+ for (sn = 0; sn <= 1; sn++) {
+ for (n = 0; n < mc->c[sn].ncubepoints; n++) {
+ MCubePoint *mcp = mc->c[sn].cubepoints + n;
+ nn = mcp->key;
+ if (nn == -1)
+ break;
+ iix = nn % 3;
+ iz = nn / 3 % mc->nz;
+ iy = nn / (3 * mc->nz) % mc->ny;
+ ix = nn / (3 * mc->nz * mc->ny);
+ fprintf(fp, "%c%d:[%d,%d,%d,%d] (%g,[%g,%g,%g],[%g,%g,%g])\n", (sn == 0 ? 'p' : 'P'),
+ n, ix, iy, iz, iix,
+ mcp->d, mcp->pos[0], mcp->pos[1], mcp->pos[2], mcp->grad[0], mcp->grad[1], mcp->grad[2]);
+ }
+ for (n = 0; n < mc->c[sn].ntriangles; n += 3) {
+ if (mc->c[sn].triangles[n] < 0)
+ break;
+ fprintf(fp, "%c%d:(%d,%d,%d)\n", (sn == 0 ? 't' : 'T'), n / 3,
+ mc->c[sn].triangles[n], mc->c[sn].triangles[n + 1], mc->c[sn].triangles[n + 2]);
+ }
+ }
+ fclose(fp);
+ }
+
+ return retval;
+}
+
+void
+MoleculeDeallocateMCube(MCube *mcube)
+{
+ free(mcube->dp);
+ free(mcube->radii);
+ free(mcube->c[0].cubepoints);
+ free(mcube->c[0].triangles);
+ free(mcube->c[1].cubepoints);
+ free(mcube->c[1].triangles);
+ free(mcube);
+}
OSDN Git Service
