Molecule *
MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
{
+ int i, n;
MoleculeFlushFrames(mp);
MoleculeInitWithAtoms(mp2, mp->atoms, mp->natoms);
if (mp->nbonds > 0) {
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->bset = NULL;
}
if (mp->mcube != NULL) {
- free(mp->mcube->dp);
- free(mp->mcube->radii);
- free(mp->mcube->c[0].cubepoints);
- free(mp->mcube->c[0].triangles);
- free(mp->mcube->c[1].cubepoints);
- free(mp->mcube->c[1].triangles);
- free(mp->mcube);
+ 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;
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, "%lf", &mview_dbuf[0]) < 1)
+ if ((i == 0 && sscanf(buf, "%lf", &dbuf[0]) < 1)
|| (i == 1 && sscanf(buf, "%lf %lf %lf",
- &mview_dbuf[1], &mview_dbuf[2], &mview_dbuf[3]) < 3)
+ &dbuf[1], &dbuf[2], &dbuf[3]) < 3)
|| (i == 2 && sscanf(buf, "%lf %lf %lf %lf",
- &mview_dbuf[4], &mview_dbuf[5], &mview_dbuf[6], &mview_dbuf[7]) < 4)) {
+ &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))
- || (strcmp(comp, "atom_resolution") == 0 && (i = 11))
- || (strcmp(comp, "bond_resolution") == 0 && (i = 12))) {
- mview_ibuf[i - 1] = atoi(valp);
- } else if ((strcmp(comp, "atom_radius") == 0 && (i = 8))
- || (strcmp(comp, "bond_radius") == 0 && (i = 9))) {
- mview_dbuf[i] = strtod(valp, NULL);
- } else if (strcmp(comp, "show_periodic_image") == 0) {
- sscanf(valp, "%d %d %d %d %d %d",
- &mview_ibuf[12], &mview_ibuf[13], &mview_ibuf[14],
- &mview_ibuf[15], &mview_ibuf[16], &mview_ibuf[17]);
+ 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)
TrackballSetRotate(mp->mview->track, mview_dbuf + 4);
}
}
+*/
return 0;
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];
}
/* 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;
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;
}
-/* Get MO coefficients for idx-th MO */
+/* 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. */
bset = mol->bset;
if (bset == NULL || bset->ncomps <= 0)
return -2; /* No basis set info */
- if (idx < 0 || idx >= bset->nmos)
+ 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) {
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. */
+/* 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
-MoleculeAllocateBasisSetRecord(Molecule *mol, Int rflag, Int ne_alpha, Int ne_beta)
+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);
return 3; /* Unexpected EOF */
}
+// Normalization constant for one gaussian component
+// 1/sqrt(Integral((Y(lm)*(r^n)*exp(-a*r*r))^2, for all r = (x, y, z)))
+// where Y(lm) is a spherical harmonic function, r^n is an "additional exponent"
+// required in expanded Molden file generated by JANPA, and a is the exponent
+// of the gaussian component.
+// The function Y(lm) is assumed so that its norm equals sqrt(4*pi/(2l+1))
+// for each m in [-l..l].
+static double
+sGaussianNormalizationConstant(int l, double a, int n)
+{
+ return 1.0/(sqrt(4 * PI / (2 * l + 1.0)) * sqrt(tgamma(l + n + 1.5) / (2.0 * pow(2.0 * a, l + n + 1.5))));
+}
+
static int
sSetupGaussianCoefficients(BasisSet *bset)
{
ShellInfo *sp;
PrimInfo *pp;
- int i, j, k;
+ int i, j, k, n;
Double *dp, d;
/* Cache the contraction coefficients for efficient calculation */
dp = bset->cns;
for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
for (j = 0, pp = bset->priminfos + sp->p_idx; j < sp->nprim; j++, pp++) {
+ n = sp->add_exp;
switch (sp->sym) {
case kGTOType_S:
+ // GNC(0,a,n) * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(0, pp->A, n);
+ *dp++ = d;
+ //{ printf("type_S: %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); }
// (8 alpha^3/pi^3)^0.25 exp(-alpha r^2)
- *dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547;
+ //*dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547;
break;
case kGTOType_P:
- // (128 alpha^5/pi^3)^0.25 [x|y|z]exp(-alpha r^2)
- d = pp->C * pow(pp->A, 1.25) * 1.425410941;
+ // GNC(1,a,n) * [x|y|z] * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(1, pp->A, n);
+ //{ printf("type_P: %g %g\n", d, pp->C * pow(pp->A, 1.25) * 1.425410941); }
+ // (128 alpha^5/pi^3)^0.25 [x|y|z]exp(-alpha r^2)
+ // d = pp->C * pow(pp->A, 1.25) * 1.425410941;
*dp++ = d;
*dp++ = d;
*dp++ = d;
break;
case kGTOType_SP:
- *dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547;
- d = pp->Csp * pow(pp->A, 1.25) * 1.425410941;
+ // GNC(0,a,n) * r^n * exp(-a*r^2)
+ *dp++ = d = pp->C * sGaussianNormalizationConstant(0, pp->A, n);
+ //{ printf("type_SP(s): %g %g\n", d, pp->C * pow(pp->A, 0.75) * 0.71270547); }
+ // GNC(1,a,n) * [x|y|z] * r^n * exp(-a*r^2)
+ d = pp->Csp * sGaussianNormalizationConstant(1, pp->A, n);
+ //{ printf("type_SP(p): %g %g\n", d, pp->Csp * pow(pp->A, 1.25) * 1.425410941); }
+ //*dp++ = pp->C * pow(pp->A, 0.75) * 0.71270547;
+ //d = pp->Csp * pow(pp->A, 1.25) * 1.425410941;
*dp++ = d;
*dp++ = d;
*dp++ = d;
break;
case kGTOType_D:
+ // GNC(2,a,n) * [xx|yy|zz] * r^n * exp(-a*r^2)
+ // GNC(2,a,n) * sqrt(3) * [xy|yz|zx] * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(2, pp->A, n);
+ //{ printf("type_D[0-2]: %g %g\n", d, pp->C * pow(pp->A, 1.75) * 1.645922781); }
+ //{ printf("type_D[3-5]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); }
+ dp[0] = dp[1] = dp[2] = d;
+ dp[3] = dp[4] = dp[5] = d * sqrt(3);
// xx|yy|zz: (2048 alpha^7/9pi^3)^0.25 [xx|yy|zz]exp(-alpha r^2)
// xy|yz|zx: (2048 alpha^7/pi^3)^0.25 [xy|xz|yz]exp(-alpha r^2)
- d = pp->C * pow(pp->A, 1.75);
- dp[0] = dp[1] = dp[2] = d * 1.645922781;
- dp[3] = dp[4] = dp[5] = d * 2.850821881;
+ // d = pp->C * pow(pp->A, 1.75);
+ //dp[0] = dp[1] = dp[2] = d * 1.645922781;
+ //dp[3] = dp[4] = dp[5] = d * 2.850821881;
dp += 6;
break;
case kGTOType_D5:
+ // D(0): GNC(2,a,n) * (1/2) * (3zz-rr) * r^n * exp(-a*r^2)
+ // D(+1): GNC(2,a,n) * sqrt(3) * xz * r^n * exp(-a*r^2)
+ // D(-1): GNC(2,a,n) * sqrt(3) * yz * r^n * exp(-a*r^2)
+ // D(+2): GNC(2,a,n) * (sqrt(3)/2) * (xx-yy) * r^n * exp(-a*r^2)
+ // D(-2): GNC(2,a,n) * sqrt(3) * xy * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(2, pp->A, n);
+ //{ printf("type_D5[0]: %g %g\n", d * 0.5, pp->C * pow(pp->A, 1.75) * 0.822961390); }
+ //{ printf("type_D5[1,2,4]: %g %g\n", d * sqrt(3), pp->C * pow(pp->A, 1.75) * 2.850821881); }
+ //{ printf("type_D5[3]: %g %g\n", d * sqrt(3) * 0.5, pp->C * pow(pp->A, 1.75) * 1.425410941); }
+ dp[0] = d * 0.5;
+ dp[1] = dp[2] = dp[4] = d * sqrt(3);
+ dp[3] = d * sqrt(3) * 0.5;
// 3zz-rr: (128 alpha^7/9pi^3)^0.25 (3zz-rr)exp(-alpha r^2)
// xy|yz|zx: (2048 alpha^7/pi^3)^0.25 [xy|xz|yz]exp(-alpha r^2)
// xx-yy: (128 alpha^7/pi^3)^0.25 (xx-yy)exp(-alpha r^2)
- d = pp->C * pow(pp->A, 1.75);
- dp[0] = d * 0.822961390;
- dp[1] = dp[2] = dp[4] = d * 2.850821881;
- dp[3] = d * 1.425410941;
+ //d = pp->C * pow(pp->A, 1.75);
+ //dp[0] = d * 0.822961390;
+ //dp[1] = dp[2] = dp[4] = d * 2.850821881;
+ //dp[3] = d * 1.425410941;
dp += 5;
break;
- /* TODO: Support F/F7 and G/G9 type orbitals */
+ case kGTOType_F:
+ // GNC(3,a,n) * [xxx|yyy|zzz] * r^n * exp(-a*r^2)
+ // GNC(3,a,n) * sqrt(5) * [xyy|xxy|xxz|xzz|yzz|yyz] * r^n * exp(-a*r^2)
+ // GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(3, pp->A, n);
+ dp[0] = dp[1] = dp[2] = d;
+ dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(5);
+ dp[9] = d * sqrt(15);
+ dp += 10;
+ break;
+ case kGTOType_F7:
+ // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2)
+ // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2)
+ // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2)
+ // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2)
+ // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
+ // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2)
+ // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(3, pp->A, n);
+ dp[0] = d * 0.5;
+ dp[1] = dp[2] = d * sqrt(3/8.0);
+ dp[3] = d * sqrt(15/4.0);
+ dp[4] = d * sqrt(15);
+ dp[5] = dp[6] = d * sqrt(5/8.0);
+ dp += 7;
+ break;
+ case kGTOType_G:
+ // GNC(4,a,n) * [xxxx|yyyy|zzzz] * exp(-a*r^2)
+ // GNC(4,a,n) * sqrt(7) * [xxxy|xxxz|yyyx|yyyz|zzzx|zzzy] * exp(-a*r^2)
+ // GNC(4,a,n) * sqrt(35/3) * [xxyy|xxzz|yyzz] * exp(-a*r^2)
+ // GNC(4,a,n) * sqrt(35) * [xxyz|yyzx|zzxy] * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(4, pp->A, n);
+ dp[0] = dp[1] = dp[2] = d;
+ dp[3] = dp[4] = dp[5] = dp[6] = dp[7] = dp[8] = d * sqrt(7);
+ dp[9] = dp[10] = dp[11] = d * sqrt(35/3.0);
+ dp[12] = dp[13] = dp[14] = d * sqrt(35);
+ dp += 15;
+ break;
+ case kGTOType_G9:
+ // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * exp(-a*r^2)
+ // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * exp(-a*r^2)
+ // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * exp(-a*r^2)
+ // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * exp(-a*r^2)
+ // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * exp(-a*r^2)
+ // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * exp(-a*r^2)
+ // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * exp(-a*r^2)
+ // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * exp(-a*r^2)
+ // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * exp(-a*r^2)
+ d = pp->C * sGaussianNormalizationConstant(4, pp->A, n);
+ dp[0] = d * 0.125;
+ dp[1] = dp[2] = d * sqrt(5/8.0);
+ dp[3] = d * sqrt(5/16.0);
+ dp[4] = d * sqrt(5/4.0);
+ dp[5] = dp[6] = d * sqrt(35/8.0);
+ dp[7] = d * sqrt(35/64.0);
+ dp[8] = d * sqrt(35/4.0);
+ dp += 9;
+ break;
}
}
}
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);
FILE *fp;
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, "\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;
}
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);
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 position. */
/* If limit is negative, its absolute value denotes the threshold distance in angstrom; otherwise,
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,
}
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;
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;
}
+int
+MoleculeReorderFrames(Molecule *mp, const Int *old_idx)
+{
+ 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;
+ }
+ 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];
+ }
+ }
+ 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];
+ }
+ }
+ 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++;
+ }
+ IntGroupIteratorRelease(&iter);
+ return n;
+}
+
#pragma mark ====== Pi Atoms ======
static inline void
#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(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
Double val, tval, *cnp, *tmpp, *mobasep, *mop;
Int i, j;
/* Cache dr and |dr|^2 */
+ if (index == 0)
+ index = bset->nmos + 1;
for (i = 0; i < mp->natoms; i++) {
Vector r;
r = ATOM_AT_INDEX(mp->atoms, i)->r;
val = 0.0;
mobasep = bset->mo + (index - 1) * bset->ncomps;
for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
- pp = bset->priminfos + sp->p_idx;
+ Double rn;
+ 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;
+ if (sp->add_exp == 0)
+ rn = 1.0;
+ else
+ rn = pow(tmpp[3], sp->add_exp * 0.5);
switch (sp->sym) {
case kGTOType_S: {
tval = 0;
for (j = 0; j < sp->nprim; j++) {
- tval += *cnp++ * exp(-pp->A * tmpp[3]);
+ tval += *cnp++ * rn * exp(-pp->A * tmpp[3]);
pp++;
}
val += mop[0] * tval;
Double x, y, z;
x = y = z = 0;
for (j = 0; j < sp->nprim; j++) {
- tval = exp(-pp->A * tmpp[3]);
+ tval = rn * exp(-pp->A * tmpp[3]);
x += *cnp++ * tval;
y += *cnp++ * tval;
z += *cnp++ * tval;
Double t, x, y, z;
t = x = y = z = 0;
for (j = 0; j < sp->nprim; j++) {
- tval = exp(-pp->A * tmpp[3]);
+ tval = rn * exp(-pp->A * tmpp[3]);
t += *cnp++ * tval;
x += *cnp++ * tval;
y += *cnp++ * tval;
Double xx, yy, zz, xy, xz, yz;
xx = yy = zz = xy = xz = yz = 0;
for (j = 0; j < sp->nprim; j++) {
- tval = exp(-pp->A * tmpp[3]);
+ tval = rn * exp(-pp->A * tmpp[3]);
xx += *cnp++ * tval;
yy += *cnp++ * tval;
zz += *cnp++ * tval;
Double d0, d1p, d1n, d2p, d2n;
d0 = d1p = d1n = d2p = d2n = 0;
for (j = 0; j < sp->nprim; j++) {
- tval = exp(-pp->A * tmpp[3]);
+ tval = rn * exp(-pp->A * tmpp[3]);
d0 += *cnp++ * tval;
d1p += *cnp++ * tval;
d1n += *cnp++ * tval;
val += d0 + d1p + d1n + d2p + d2n;
break;
}
- /* TODO: Support F/F7 and G/G9 type orbitals */
+ case kGTOType_F: {
+ Double xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz;
+ xxx = yyy = zzz = xyy = xxy = xxz = xzz = yzz = yyz = xyz = 0;
+ for (j = 0; j < sp->nprim; j++) {
+ tval = rn * exp(-pp->A * tmpp[3]);
+ xxx += *cnp++ * tval;
+ yyy += *cnp++ * tval;
+ zzz += *cnp++ * tval;
+ xyy += *cnp++ * tval;
+ xxy += *cnp++ * tval;
+ xxz += *cnp++ * tval;
+ xzz += *cnp++ * tval;
+ yzz += *cnp++ * tval;
+ yyz += *cnp++ * tval;
+ xyz += *cnp++ * tval;
+ pp++;
+ }
+ xxx *= mop[0] * tmpp[0] * tmpp[0] * tmpp[0];
+ yyy *= mop[1] * tmpp[1] * tmpp[1] * tmpp[1];
+ zzz *= mop[2] * tmpp[2] * tmpp[2] * tmpp[2];
+ xyy *= mop[3] * tmpp[0] * tmpp[1] * tmpp[1];
+ xxy *= mop[4] * tmpp[0] * tmpp[0] * tmpp[1];
+ xxz *= mop[5] * tmpp[0] * tmpp[0] * tmpp[2];
+ xzz *= mop[6] * tmpp[0] * tmpp[2] * tmpp[2];
+ yzz *= mop[7] * tmpp[1] * tmpp[2] * tmpp[2];
+ yyz *= mop[8] * tmpp[1] * tmpp[1] * tmpp[2];
+ xyz *= mop[9] * tmpp[0] * tmpp[1] * tmpp[2];
+ val += xxx + yyy + zzz + xyy + xxy + xxz + xzz + yzz + yyz + xyz;
+ break;
+ }
+ case kGTOType_F7: {
+ Double f0, f1p, f1n, f2p, f2n, f3p, f3n;
+ f0 = f1p = f1n = f2p = f2n = f3p = f3n = 0;
+ for (j = 0; j < sp->nprim; j++) {
+ tval = rn * exp(-pp->A * tmpp[3]);
+ f0 += *cnp++ * tval;
+ f1p += *cnp++ * tval;
+ f1n += *cnp++ * tval;
+ f2p += *cnp++ * tval;
+ f2n += *cnp++ * tval;
+ f3p += *cnp++ * tval;
+ f3n += *cnp++ * tval;
+ pp++;
+ }
+ // F(0): GNC(3,a,n) * (1/2) * (5zzz-3zrr) * r^n * exp(-a*r^2)
+ // F(+1): GNC(3,a,n) * sqrt(3/8) * (5xzz-xrr) * r^n * exp(-a*r^2)
+ // F(-1): GNC(3,a,n) * sqrt(3/8) * (5yzz-yrr) * r^n * exp(-a*r^2)
+ // F(+2): GNC(3,a,n) * sqrt(15/4) * (xxz-yyz) * r^n * exp(-a*r^2)
+ // F(-2): GNC(3,a,n) * sqrt(15) * xyz * r^n * exp(-a*r^2)
+ // F(+3): GNC(3,a,n) * sqrt(5/8) * (xxx-3xyy) * r^n * exp(-a*r^2)
+ // F(-3): GNC(3,a,n) * sqrt(5/8) * (3xxy-yyy) * r^n * exp(-a*r^2)
+ f0 *= mop[0] * tmpp[2] * (5 * tmpp[2] * tmpp[2] - 3 * tmpp[3]);
+ f1p *= mop[1] * tmpp[0] * (5 * tmpp[2] * tmpp[2] - tmpp[3]);
+ f1n *= mop[2] * tmpp[1] * (5 * tmpp[2] * tmpp[2] - tmpp[3]);
+ f2p *= mop[3] * tmpp[2] * (tmpp[0] * tmpp[0] - tmpp[1] * tmpp[1]);
+ f2n *= mop[4] * tmpp[0] * tmpp[1] * tmpp[2];
+ f3p *= mop[5] * tmpp[0] * (tmpp[0] * tmpp[0] - 3 * tmpp[1] * tmpp[1]);
+ f3n *= mop[6] * tmpp[1] * (3 * tmpp[0] * tmpp[0] - tmpp[2] * tmpp[2]);
+ val += f0 + f1p + f1n + f2p + f2n + f3p + f3n;
+ break;
+ }
+ case kGTOType_G: {
+ Double xxxx, yyyy, zzzz, xxxy, xxxz, yyyx, yyyz, zzzx, zzzy, xxyy, xxzz, yyzz, xxyz, yyxz, zzxy;
+ xxxx = yyyy = zzzz = xxxy = xxxz = yyyx = yyyz = zzzx = zzzy = xxyy = xxzz = yyzz = xxyz = yyxz = zzxy = 0;
+ for (j = 0; j < sp->nprim; j++) {
+ tval = rn * exp(-pp->A * tmpp[3]);
+ xxxx += *cnp++ * tval;
+ yyyy += *cnp++ * tval;
+ zzzz += *cnp++ * tval;
+ xxxy += *cnp++ * tval;
+ xxxz += *cnp++ * tval;
+ yyyx += *cnp++ * tval;
+ yyyz += *cnp++ * tval;
+ zzzx += *cnp++ * tval;
+ zzzy += *cnp++ * tval;
+ xxyy += *cnp++ * tval;
+ xxzz += *cnp++ * tval;
+ yyzz += *cnp++ * tval;
+ xxyz += *cnp++ * tval;
+ yyxz += *cnp++ * tval;
+ zzxy += *cnp++ * tval;
+ pp++;
+ }
+ xxxx *= mop[0] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[0];
+ yyyy *= mop[1] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[1];
+ zzzz *= mop[2] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[2];
+ xxxy *= mop[3] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[1];
+ xxxz *= mop[4] * tmpp[0] * tmpp[0] * tmpp[0] * tmpp[2];
+ yyyx *= mop[5] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[0];
+ yyyz *= mop[6] * tmpp[1] * tmpp[1] * tmpp[1] * tmpp[2];
+ zzzx *= mop[7] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[0];
+ zzzy *= mop[8] * tmpp[2] * tmpp[2] * tmpp[2] * tmpp[1];
+ xxyy *= mop[9] * tmpp[0] * tmpp[0] * tmpp[1] * tmpp[1];
+ xxzz *= mop[10] * tmpp[0] * tmpp[0] * tmpp[2] * tmpp[2];
+ yyzz *= mop[11] * tmpp[1] * tmpp[1] * tmpp[2] * tmpp[2];
+ xxyz *= mop[12] * tmpp[0] * tmpp[0] * tmpp[1] * tmpp[2];
+ yyxz *= mop[13] * tmpp[1] * tmpp[1] * tmpp[0] * tmpp[2];
+ zzxy *= mop[14] * tmpp[2] * tmpp[2] * tmpp[0] * tmpp[1];
+ val += xxxx + yyyy + zzzz + xxxy + xxxz + yyyx + yyyz + zzzx + zzzy + xxyy + xxzz + yyzz + xxyz + yyxz + zzxy;
+ break;
+ }
+ case kGTOType_G9: {
+ Double g0, g1p, g1n, g2p, g2n, g3p, g3n, g4p, g4n;
+ Double xx = tmpp[0] * tmpp[0];
+ Double yy = tmpp[1] * tmpp[1];
+ Double zz = tmpp[2] * tmpp[2];
+ Double rr = tmpp[3];
+ g0 = g1p = g1n = g2p = g2n = g3p = g3n = g4p = g4n = 0;
+ for (j = 0; j < sp->nprim; j++) {
+ tval = rn * exp(-pp->A * tmpp[3]);
+ g0 += *cnp++ * tval;
+ g1p += *cnp++ * tval;
+ g1n += *cnp++ * tval;
+ g2p += *cnp++ * tval;
+ g2n += *cnp++ * tval;
+ g3p += *cnp++ * tval;
+ g3n += *cnp++ * tval;
+ g4p += *cnp++ * tval;
+ g4n += *cnp++ * tval;
+ pp++;
+ }
+ // G(0): GNC(4,a,n) * (1/8) * (35zzzz-30zzrr+3rrrr) * r^n * exp(-a*r^2)
+ // G(+1): GNC(4,a,n) * sqrt(5/8) * (7xzzz-3xzrr) * r^n * exp(-a*r^2)
+ // G(-1): GNC(4,a,n) * sqrt(5/8) * (7yzzz-3yzrr) * r^n * exp(-a*r^2)
+ // G(+2): GNC(4,a,n) * sqrt(5/16) * (xx-yy)(7zz-rr) * r^n * exp(-a*r^2)
+ // G(-2): GNC(4,a,n) * sqrt(5/4) * (7xyzz-xyrr) * r^n * exp(-a*r^2)
+ // G(+3): GNC(4,a,n) * sqrt(35/8) * (xxxz-3xyyz) * r^n * exp(-a*r^2)
+ // G(-3): GNC(4,a,n) * sqrt(35/8) * (3xxyz-yyyz) * r^n * exp(-a*r^2)
+ // G(+4): GNC(4,a,n) * sqrt(35/64) * (xxxx-6xxyy+yyyy) * r^n * exp(-a*r^2)
+ // G(-4): GNC(4,a,n) * sqrt(35/4) * (xxxy-xyyy) * r^n * exp(-a*r^2)
+ g0 *= mop[0] * (35 * zz * zz - 30 * zz * rr + 3 * rr * rr);
+ g1p *= mop[1] * tmpp[0] * tmpp[2] * (7 * zz - 3 * rr);
+ g1n *= mop[2] * tmpp[1] * tmpp[2] * (7 * zz - 3 * rr);
+ g2p *= mop[3] * (xx - yy) * (7 * zz - rr);
+ g2n *= mop[4] * tmpp[0] * tmpp[1] * (7 * zz - rr);
+ g3p *= mop[5] * tmpp[0] * tmpp[2] * (xx - 3 * yy);
+ g3n *= mop[6] * tmpp[1] * tmpp[2] * (3 * xx - yy);
+ g4p *= mop[7] * (xx * xx - 6 * xx * yy + yy * yy);
+ g4n *= mop[8] * tmpp[0] * tmpp[1] * (xx - yy);
+ val += g0 + g1p + g1n + g2p + g2n + g3p + g3n + g4p + g4n;
+ break;
+ }
}
}
return val;
}
/* Calculate one MO. The input vectors are angstrom unit (changed from bohr unit: 20140520) */
-/* mono is the MO number (1-based) */
+/* 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)
{
/* 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 = (int)floor(log10(fabs(d)));
- Double base = d * pow(10, -1.0 * exponent);
+ 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)
}
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;
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;
}
int
MoleculeUpdateMCube(Molecule *mol, int idn)
{
- Int flags, retval, step, hstep, sn;
+ Int retval, step, sn;
Int n, ix, iy, iz, nx, ny, nz;
Int nn, iix, iiy, iiz;
Int ncubepoints, c1, c2, c3;
return -1; /* Number of atoms is smaller than expected */
mc = mol->mcube;
- if (idn > 0) {
+ if (idn >= 0) {
ShellInfo *sp;
Double *mobasep, *mop, mopmax;
Double xmin, xmax, ymin, ymax, zmin, zmax;
if (mc->radii == NULL)
return -2; /* Out of memory */
mc->nradii = mol->natoms;
- memset(mc->radii, 0, sizeof(Double) * mc->nradii);
- mobasep = mol->bset->mo + (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;
+ 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;
n = (ix * ny + iy) * nz + iz;
if (mc->dp[n] == DBL_MAX) {
p.z = mc->origin.z + mc->dz * iz;
- mc->dp[n] = sCalcMOPoint(mol, mol->bset, mc->idn, &p, tmp);
+ 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);
+ }
}
}
}
/* Less than 3 points: no triangles */
if (mc->c[sn].ntriangles > 0)
mc->c[sn].triangles[0] = -1; /* End mark */
- retval = 0;
- goto end;
+ continue;
}
/* Create triangle table */
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
