X-Git-Url: http://git.osdn.net/view?a=blobdiff_plain;f=MolLib%2FMolecule.c;h=ac8f541608afff5a6bd2d8a978803f4299757ac3;hb=3d3ec3d7aa41664271bdfebe8fb8566f77e1e818;hp=66c0bcfb7d722fdf1516c596ebbabe97b3c39b88;hpb=f4f5b2c186b43f6b03df2fae7cea0b4764855953;p=molby%2FMolby.git diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c index 66c0bcf..ac8f541 100755 --- a/MolLib/Molecule.c +++ b/MolLib/Molecule.c @@ -19,6 +19,7 @@ #include #include #include +#include #include "Missing.h" #include "Dcd.h" @@ -154,8 +155,8 @@ BasisSetRelease(BasisSet *bset) 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) { @@ -343,6 +344,7 @@ MoleculeInitWithAtoms(Molecule *mp, const Atom *atoms, int natoms) Molecule * MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp) { + int i, n; MoleculeFlushFrames(mp); MoleculeInitWithAtoms(mp2, mp->atoms, mp->natoms); if (mp->nbonds > 0) { @@ -386,6 +388,17 @@ MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp) 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) @@ -467,6 +480,7 @@ MoleculeRetain(Molecule *mp) void MoleculeClear(Molecule *mp) { + int i; if (mp == NULL) return; if (mp->arena != NULL) { @@ -477,12 +491,7 @@ MoleculeClear(Molecule *mp) 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); @@ -536,6 +545,19 @@ MoleculeClear(Molecule *mp) 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; @@ -766,8 +788,8 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -788,9 +810,9 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -1563,14 +1585,13 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -1581,15 +1602,23 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -1599,6 +1628,8 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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) { @@ -1606,24 +1637,334 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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 */ @@ -1634,7 +1975,8 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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) @@ -1655,19 +1997,28 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -1725,17 +2076,6 @@ MoleculeLoadPsfFile(Molecule *mp, const char *fname, char **errbuf) 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++) { @@ -1945,10 +2285,9 @@ MoleculeLoadPsfFile(Molecule *mp, const char *fname, char **errbuf) 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]; } @@ -2159,7 +2498,7 @@ MoleculeLoadTepFile(Molecule *mp, const char *fname, char **errbuf) } } 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); @@ -2357,7 +2696,7 @@ MoleculeLoadShelxFile(Molecule *mp, const char *fname, char **errbuf) 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); @@ -2371,7 +2710,7 @@ MoleculeLoadShelxFile(Molecule *mp, const char *fname, char **errbuf) /* 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; @@ -2408,6 +2747,10 @@ MoleculeAddGaussianOrbitalShell(Molecule *mol, Int sym, Int nprims, Int a_idx) 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; } @@ -2434,7 +2777,86 @@ MoleculeAddGaussianPrimitiveCoefficients(Molecule *mol, Double exponent, Double 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) { @@ -2463,8 +2885,8 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou 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); @@ -2476,8 +2898,14 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou 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) @@ -2492,36 +2920,63 @@ MoleculeSetMOCoefficients(Molecule *mol, Int idx, Double energy, Int ncomps, Dou 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; @@ -2577,12 +3032,25 @@ sReadNumberArray(void *basep, Int *countp, Int size, Int num, FILE *fp, int *lnp 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 */ @@ -2598,44 +3066,134 @@ sSetupGaussianCoefficients(BasisSet *bset) 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; } } } @@ -2656,7 +3214,7 @@ MoleculeLoadGaussianFchkFile(Molecule *mp, const char *fname, char **errbuf) Int *iary; Double *dary; Atom *ap; - Vector *vp; +/* Vector *vp; */ Double w; *errbuf = NULL; @@ -2702,10 +3260,11 @@ MoleculeLoadGaussianFchkFile(Molecule *mp, const char *fname, char **errbuf) 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) { @@ -2788,13 +3347,10 @@ MoleculeLoadGaussianFchkFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -3212,7 +3768,7 @@ MoleculeLoadGamessDatFile(Molecule *mol, const char *fname, char **errbuf) 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; @@ -3247,7 +3803,7 @@ MoleculeLoadGamessDatFile(Molecule *mol, const char *fname, char **errbuf) } 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); @@ -3292,7 +3848,7 @@ exit_loop: 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); @@ -3816,6 +4372,7 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; @@ -4137,7 +4694,7 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf) } 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"); @@ -4164,9 +4721,141 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf) 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; } @@ -4275,24 +4964,6 @@ MoleculeWriteToPsfFile(Molecule *mp, const char *fname, char **errbuf) 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); @@ -4744,58 +5415,6 @@ sOutputBondInstructions(FILE *fp, int natoms, Atom *atoms, int overlap_correctio 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) @@ -5126,7 +5745,10 @@ MoleculeDeserialize(const char *data, Int length, Int *timep) 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); @@ -5685,25 +6307,26 @@ MoleculeSearchImpropersAcrossAtomGroup(Molecule *mp, IntGroup *atomgroup) 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) @@ -5713,11 +6336,11 @@ MoleculeFindCloseAtoms(Molecule *mp, Int index, Double limit, Int *outNbonds, In 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; @@ -5735,42 +6358,19 @@ int 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; @@ -6543,9 +7143,10 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice ap = ATOM_AT_INDEX(mp->atoms, n); if (SYMOP_ALIVE(ap->symop)) { /* Calculate the cumulative symop */ + Transform tr2; MoleculeGetTransformForSymop(mp, ap->symop, &t1, 0); - TransformMul(tr, tr, t1); - if (MoleculeGetSymopForTransform(mp, tr, &symop1, 0) != 0) { + TransformMul(tr2, tr, t1); + if (MoleculeGetSymopForTransform(mp, tr2, &symop1, 0) != 0) { if (indices != NULL) indices[i] = -1; continue; /* Skip this atom */ @@ -6584,9 +7185,9 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice 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) @@ -6600,20 +7201,22 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice for (i = n0; i < n1; i++) { Int b[2], j; ap = ATOM_AT_INDEX(mp->atoms, i); - if (MoleculeGetSymopForTransform(mp, tr, &ap->symop, 0) == 0) { + 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++) { - nr = ATOM_AT_INDEX(mp->atoms, cp[n])->r; + 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 == ap->atomicNumber && VecLength2(dr) < 1e-6) + if (ap2->atomicNumber == apn->atomicNumber && VecLength2(dr) < 1e-6) break; } if (j < mp->natoms) { @@ -6626,28 +7229,6 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice } } } - -#if 0 - for (i = 0; i < n0; i++) { - int b[2]; - Int *cp; - b[0] = table[i]; - if (b[0] < 0 || b[0] == i) - continue; - 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; - if (MoleculeLookupBond(mp, b[0], b[1]) >= 0) - continue; - MoleculeAddBonds(mp, 1, b, NULL, 1); - } - } -#endif mp->needsMDRebuild = 1; __MoleculeUnlock(mp); free(table); @@ -6870,7 +7451,9 @@ int 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 */ @@ -8359,7 +8942,7 @@ MoleculeDeleteAngles(Molecule *mp, Int *angles, IntGroup *where) __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) { @@ -8411,7 +8994,7 @@ MoleculeDeleteDihedrals(Molecule *mp, Int *dihedrals, IntGroup *where) __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) { @@ -8463,7 +9046,7 @@ MoleculeDeleteImpropers(Molecule *mp, Int *impropers, IntGroup *where) __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) { @@ -8908,12 +9491,12 @@ MoleculeChangeResidueNumberWithArray(Molecule *mp, IntGroup *group, Int *resSeqs 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); @@ -8954,7 +9537,7 @@ MoleculeChangeResidueNumberWithArray(Molecule *mp, IntGroup *group, Int *resSeqs 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, @@ -9805,6 +10388,7 @@ MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double } 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; @@ -10134,6 +10718,9 @@ MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const 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; @@ -10157,17 +10744,14 @@ MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const /* 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; @@ -10181,6 +10765,18 @@ MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const } } + /* 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, @@ -10190,23 +10786,28 @@ MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const ... 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); @@ -10223,23 +10824,37 @@ MoleculeInsertFrames(Molecule *mp, IntGroup *group, const Vector *inFrame, const 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; @@ -10255,6 +10870,7 @@ MoleculeRemoveFrames(Molecule *mp, IntGroup *inGroup, Vector *outFrame, Vector * 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) @@ -10320,38 +10936,43 @@ MoleculeRemoveFrames(Molecule *mp, IntGroup *inGroup, Vector *outFrame, Vector * 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); @@ -10360,18 +10981,20 @@ MoleculeRemoveFrames(Molecule *mp, IntGroup *inGroup, Vector *outFrame, Vector * 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) @@ -10383,19 +11006,29 @@ MoleculeRemoveFrames(Molecule *mp, IntGroup *inGroup, Vector *outFrame, Vector * 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; @@ -10429,7 +11062,7 @@ MoleculeSelectFrame(Molecule *mp, int frame, int copyback) /* 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; @@ -10446,7 +11079,7 @@ MoleculeSelectFrame(Molecule *mp, int frame, int copyback) 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); @@ -10471,25 +11104,205 @@ MoleculeFlushFrames(Molecule *mp) 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 +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 +MoleculeUpdatePiAnchorPositions(Molecule *mol) +{ + 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, *ap2; - Int i, n, *ip; + Atom *ap; if (mol == NULL || idx < 0 || idx >= mol->natoms) return; ap = ATOM_AT_INDEX(mol->atoms, idx); if (ap->anchor == NULL) return; - ip = AtomConnectData(&ap->anchor->connect); - n = ap->anchor->connect.count; - VecZero(ap->r); - for (i = 0; i < ap->anchor->connect.count; i++) { - ap2 = ATOM_AT_INDEX(mol->atoms, ip[i]); - VecScaleInc(ap->r, ap2->r, ap->anchor->coeffs[i]); - } + sMoleculeCalculatePiAnchorPosition(ap, mol->atoms); } int @@ -10626,36 +11439,46 @@ MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Doub #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 */ 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; @@ -10665,7 +11488,7 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp) 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; @@ -10681,7 +11504,7 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp) 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; @@ -10699,7 +11522,7 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp) 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; @@ -10721,7 +11544,7 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp) 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; @@ -10737,14 +11560,155 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp) 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 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) { @@ -10757,6 +11721,9 @@ MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, cons 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; @@ -10776,7 +11743,7 @@ MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, cons 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); @@ -10790,7 +11757,7 @@ MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, cons 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; @@ -10809,35 +11776,38 @@ MoleculeCalcMO(Molecule *mp, Int mono, const Vector *op, const Vector *dxp, cons 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; @@ -10914,17 +11884,21 @@ MoleculeOutputCube(Molecule *mp, Int index, const char *fname, const char *comme 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); @@ -10932,7 +11906,20 @@ MoleculeOutputCube(Molecule *mp, Int index, const char *fname, const char *comme 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"); } @@ -10941,3 +11928,935 @@ MoleculeOutputCube(Molecule *mp, Int index, const char *fname, const char *comme 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); +}