X-Git-Url: http://git.osdn.net/view?a=blobdiff_plain;f=MolLib%2FMolecule.c;h=4672586fb56c22eb4579281e7fe88a7e9118f38f;hb=388f1f455139811b1c290d55131fcc1aebf6a5cd;hp=118991db63e32e0225ae8abfbd2b9279e75d3f2a;hpb=cb32e5f43a082f714ffc00e8f6a07b8638937390;p=molby%2FMolby.git

diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 118991d..4672586 100755
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -19,6 +19,7 @@
 #include <stddef.h>
 #include <ctype.h>
 #include <math.h>
+#include <float.h>
 
 #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) {
@@ -258,34 +259,6 @@ AtomConnectHasEntry(AtomConnect *ac, Int ent)
 	return 0;
 }
 
-#if PIATOM
-void
-PiAtomDuplicate(PiAtom *pa, const PiAtom *cpa)
-{
-	memmove(pa, cpa, sizeof(PiAtom));
-	pa->connect.count = 0;
-	AtomConnectResize(&(pa->connect), cpa->connect.count);
-	memmove(AtomConnectData(&(pa->connect)), AtomConnectData((AtomConnect *)&(cpa->connect)), sizeof(Int) * cpa->connect.count);
-	pa->ncoeffs = 0;
-	pa->coeffs = NULL;
-	if (cpa->ncoeffs > 0) {
-		NewArray(&(pa->coeffs), &(pa->ncoeffs), sizeof(Double), cpa->ncoeffs);
-		memmove(pa->coeffs, cpa->coeffs, sizeof(Double) * cpa->ncoeffs);
-	}
-}
-
-void
-PiAtomClean(PiAtom *pa)
-{
-	AtomConnectResize(&(pa->connect), 0);
-	pa->ncoeffs = 0;
-	if (pa->coeffs != NULL) {
-		free(pa->coeffs);
-		pa->coeffs = NULL;
-	}
-}
-#endif
-
 #pragma mark ====== Accessor types ======
 
 MolEnumerable *
@@ -371,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) {
@@ -406,29 +380,25 @@ MoleculeInitWithMolecule(Molecule *mp2, Molecule *mp)
 		NewArray(&(mp2->syms), &(mp2->nsyms), sizeof(Transform), mp->nsyms);
 		memmove(mp2->syms, mp->syms, sizeof(Transform) * mp2->nsyms);
 	}
-#if PIATOM
-	if (mp->npiatoms > 0) {
-		NewArray(&(mp2->piatoms), &(mp2->npiatoms), sizeof(PiAtom), mp->npiatoms);
-		for (i = 0; i < mp->npiatoms; i++)
-			PiAtomDuplicate(mp2->piatoms + i, mp->piatoms + i);
-	}
-	if (mp->npibonds > 0) {
-		NewArray(&(mp2->pibonds), &(mp2->npibonds), sizeof(Int) * 4, mp->npibonds);
-		memmove(mp2->pibonds, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
-	}
-	if (mp->npiconnects > 0) {
-		NewArray(&(mp2->piconnects), &(mp2->npiconnects), sizeof(Int), mp->npiconnects);
-		memmove(mp2->piconnects, mp->piconnects, sizeof(Int) * mp->npiconnects);
-	}
-#endif
-	
-/*	mp2->useFlexibleCell = mp->useFlexibleCell; */
+
+	/*	mp2->useFlexibleCell = mp->useFlexibleCell; */
 	if (mp->nframe_cells > 0) {
 		if (NewArray(&mp2->frame_cells, &mp2->nframe_cells, sizeof(Vector) * 4, mp->nframe_cells) == NULL)
 			goto error;
 		memmove(mp2->frame_cells, mp->frame_cells, sizeof(Vector) * 4 * mp->nframe_cells);
 	}
 	
+	if (mp->nmolprops > 0) {
+		if (NewArray(&mp2->molprops, &mp2->nmolprops, sizeof(MolProp), mp->nmolprops) == NULL)
+			goto error;
+		n = MoleculeGetNumberOfFrames(mp);
+		for (i = 0; i < mp2->nmolprops; i++) {
+			mp2->molprops[i].propname = strdup(mp->molprops[i].propname);
+			mp2->molprops[i].propvals = (Double *)malloc(sizeof(Double) * n);
+			memcpy(mp2->molprops[i].propvals, mp->molprops[i].propvals, sizeof(Double) * n);
+		}
+	}
+	
 	/* FIXME: should bset (basis set info) and elpot be duplicated or not?  */
 
 	if (mp->par != NULL)
@@ -510,6 +480,7 @@ MoleculeRetain(Molecule *mp)
 void
 MoleculeClear(Molecule *mp)
 {
+	int i;
 	if (mp == NULL)
 		return;
 	if (mp->arena != NULL) {
@@ -520,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);
@@ -566,26 +532,6 @@ MoleculeClear(Molecule *mp)
 		mp->syms = NULL;
 		mp->nsyms = 0;
 	}
-#if PIATOM
-	if (mp->piatoms != NULL) {
-		for (i = 0; i < mp->npiatoms; i++) {
-			PiAtomClean(mp->piatoms + i);
-		}
-		free(mp->piatoms);
-		mp->piatoms = NULL;
-		mp->npiatoms = 0;
-	}
-	if (mp->pibonds != NULL) {
-		free(mp->pibonds);
-		mp->pibonds = NULL;
-		mp->npibonds = 0;
-	}
-	if (mp->piconnects != NULL) {
-		free(mp->piconnects);
-		mp->piconnects = NULL;
-		mp->npiconnects = 0;
-	}
-#endif
 	if (mp->selection != NULL) {
 		IntGroupRelease(mp->selection);
 		mp->selection = NULL;
@@ -599,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;
@@ -829,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;
@@ -851,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;
@@ -939,6 +898,27 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 				i++;
 			}
 			continue;
+		} else if (strcmp(buf, "!:uff_types") == 0) {
+			i = 0;
+			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+				if (buf[0] == '!')
+					continue;
+				if (buf[0] == '\n')
+					break;
+				/* idx uff_type */
+				if (sscanf(buf, "%d %6s", &ibuf[0], cbuf[0]) < 2) {
+					s_append_asprintf(errbuf, "line %d: uff type info cannot be read for atom %d", lineNumber, i + 1);
+					goto err_exit;
+				}
+				if (i >= mp->natoms) {
+					s_append_asprintf(errbuf, "line %d: too many uff type info\n", lineNumber);
+					goto err_exit;
+				}
+				ap = ATOM_AT_INDEX(mp->atoms, i);
+				strncpy(ap->uff_type, cbuf[0], 5);
+				ap->uff_type[5] = 0;
+				i++;
+			}
 		} else if (strcmp(buf, "!:mm_exclude") == 0) {
 			i = 0;
 			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
@@ -1089,6 +1069,36 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 				}
 			}
 			continue;
+		} else if (strcmp(buf, "!:bond_orders") == 0) {
+			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
+				if (buf[0] == '!')
+					continue;
+				if (buf[0] == '\n')
+					break;
+				/* b1 b2 b3 b4 */
+				i = sscanf(buf, "%lf %lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2], &dbuf[3]);
+				if (i == 0) {
+					s_append_asprintf(errbuf, "line %d: bad bond order format", lineNumber);
+					goto err_exit;
+				}
+				for (j = 0; j < i; j++) {
+					AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbondOrders, &dbuf[j]);
+				}
+			}
+			if (mp->nbondOrders > mp->nbonds) {
+				s_append_asprintf(errbuf, "line %d: warning: the number of bond order info (%d) exceeds number of bonds (%d) - ignoring excess info\n", lineNumber, mp->nbondOrders, mp->nbonds);
+				nwarnings++;
+				mp->nbondOrders = mp->nbonds;
+			} else if (mp->nbondOrders < mp->nbonds) {
+				s_append_asprintf(errbuf, "line %d: warning: the number of bond order info (%d) is less than number of bonds (%d)\n", lineNumber, mp->nbondOrders, mp->nbonds);
+				nwarnings++;
+				j = mp->nbondOrders;
+				AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbonds - 1, NULL);
+				for (i = j; i < mp->nbonds; i++)
+					mp->bondOrders[i] = 0.0;
+			}
+			continue;
+			
 		} else if (strcmp(buf, "!:angles") == 0) {
 			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
 				if (buf[0] == '!')
@@ -1108,7 +1118,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 					if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2]) {
 						s_append_asprintf(errbuf, "line %d: warning: bad angle specification (%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2]);
 						nwarnings++;
-					} else if (MoleculeAreAtomsConnected(mp, iibuf[0], iibuf[1]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0) {
+					} else if (MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0) {
 						s_append_asprintf(errbuf, "line %d: warning: angle with non-bonded atoms (%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2]);
 						nwarnings++;						
 					} else if (MoleculeLookupAngle(mp, iibuf[0], iibuf[1], iibuf[2]) >= 0) {
@@ -1140,7 +1150,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 					if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[3] < 0 || iibuf[3] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2] || iibuf[2] == iibuf[3] || iibuf[0] == iibuf[2] || iibuf[1] == iibuf[3] || iibuf[0] == iibuf[3]) {
 						s_append_asprintf(errbuf, "line %d: warning: bad dihedral specification (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
 						nwarnings++;
-					} else if (MoleculeAreAtomsConnected(mp, iibuf[0], iibuf[1]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
+					} else if (MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
 						s_append_asprintf(errbuf, "line %d: warning: dihedral with non-bonded atoms (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
 						nwarnings++;						
 					} else if (MoleculeLookupDihedral(mp, iibuf[0], iibuf[1], iibuf[2], iibuf[3]) >= 0) {
@@ -1172,7 +1182,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 					if (iibuf[0] < 0 || iibuf[0] >= mp->natoms || iibuf[1] < 0 || iibuf[1] >= mp->natoms || iibuf[2] < 0 || iibuf[2] >= mp->natoms || iibuf[3] < 0 || iibuf[3] >= mp->natoms || iibuf[0] == iibuf[1] || iibuf[1] == iibuf[2] || iibuf[2] == iibuf[3] || iibuf[0] == iibuf[2] || iibuf[1] == iibuf[3] || iibuf[0] == iibuf[3]) {
 						s_append_asprintf(errbuf, "line %d: warning: bad improper specification (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
 						nwarnings++;
-					} else if (MoleculeAreAtomsConnected(mp, iibuf[0], iibuf[2]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[1], iibuf[2]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
+					} else if (MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[0]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[1]) == 0 || MoleculeAreAtomsConnected(mp, iibuf[2], iibuf[3]) == 0) {
 						s_append_asprintf(errbuf, "line %d: warning: improper with non-bonded atoms (%d-%d-%d-%d) - skipped\n", lineNumber, iibuf[0], iibuf[1], iibuf[2], iibuf[3]);
 						nwarnings++;						
 					} else if (MoleculeLookupImproper(mp, iibuf[0], iibuf[1], iibuf[2], iibuf[3]) >= 0) {
@@ -1245,69 +1255,6 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 				}
 			}
 			continue;
-#if PIATOM
-		} else if (strcmp(buf, "!:pi_atoms") == 0) {
-			PiAtom *pp;
-			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
-				if (buf[0] == '!')
-					continue;
-				if (buf[0] == '\n')
-					break;
-				if (sscanf(buf, "%6s %6s %d", cbuf[0], cbuf[1], &ibuf[0]) < 3) {
-					s_append_asprintf(errbuf, "line %d: pi atoms info cannot be read", lineNumber);
-					goto err_exit;
-				}
-				pp = (PiAtom *)AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, NULL);
-				memset(pp, 0, sizeof(PiAtom));
-				strncpy(pp->aname, cbuf[0], 4);
-				pp->type = AtomTypeEncodeToUInt(cbuf[1]);
-				if (ibuf[0] <= 0)
-					continue;
-				AtomConnectResize(&pp->connect, ibuf[0]);
-				ip = AtomConnectData(&pp->connect);
-				NewArray(&pp->coeffs, &pp->ncoeffs, sizeof(Double), ibuf[0]);
-				for (i = 0; i < ibuf[0]; i++) {
-					if (ReadLine(buf, sizeof buf, fp, &lineNumber) <= 0) {
-						s_append_asprintf(errbuf, "line %d: unexpected end of file during reading pi atoms info", lineNumber);
-						goto err_exit;
-					}
-					if (sscanf(buf, "%d %lf", &ibuf[1], &dbuf[0]) < 2) {
-						s_append_asprintf(errbuf, "line %d: bad format during pi atoms info", lineNumber);
-						goto err_exit;
-					}
-					ip[i] = ibuf[1];
-					pp->coeffs[i] = dbuf[0];
-				}
-			}
-			continue;
-		} else if (strcmp(buf, "!:pi_atom_constructs") == 0) {
-			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
-				if (buf[0] == '!')
-					continue;
-				if (buf[0] == '\n')
-					break;
-				/* a1 b1 c1 d1 a2 b2 c2 d2 */ 
-				i = sscanf(buf, "%d %d %d %d %d %d %d %d", &ibuf[0], &ibuf[1], &ibuf[2], &ibuf[3], &ibuf[4], &ibuf[5], &ibuf[6], &ibuf[7]);
-				if (i == 0 || i % 4 != 0)
-					goto pi_atom_constructs_bad_format;
-				for (j = 0; j < i; j++) {
-					if (ibuf[j] <= -2) {
-						ibuf[j] = (-ibuf[j] - 2) + ATOMS_MAX_NUMBER;
-						if (ibuf[j] - ATOMS_MAX_NUMBER >= mp->npiatoms)
-							goto pi_atom_constructs_bad_format;
-					} else if (ibuf[j] >= mp->natoms) {
-						goto pi_atom_constructs_bad_format;
-					}
-					if (j % 4 == 3) {
-						AssignArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, mp->npibonds, &ibuf[j - 3]);
-					}
-				}
-			}
-			continue;
-		pi_atom_constructs_bad_format:
-			s_append_asprintf(errbuf, "line %d: bad format in pi_atom_constructs", lineNumber);
-			goto err_exit;
-#endif  /*  PIATOM  */
 		} else if (strcmp(buf, "!:anisotropic_thermal_parameters") == 0) {
 			i = 0;
 			while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
@@ -1518,6 +1465,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 						   || (strcmp(comp, "cutoff") == 0 && (dp = &arena->cutoff) != NULL)
 						   || (strcmp(comp, "electro_cutoff") == 0 && (dp = &arena->electro_cutoff) != NULL)
 						   || (strcmp(comp, "pairlist_distance") == 0 && (dp = &arena->pairlist_distance) != NULL)
+						   || (strcmp(comp, "switch_distance") == 0 && (dp = &arena->switch_distance) != NULL)
 						   || (strcmp(comp, "temperature") == 0 && (dp = &arena->temperature) != NULL)
 						   || (strcmp(comp, "andersen_coupling") == 0 && (dp = &arena->andersen_thermo_coupling) != NULL)
 						   || (strcmp(comp, "dielectric") == 0 && (dp = &arena->dielectric) != NULL)
@@ -1637,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;
@@ -1655,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;
@@ -1673,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) {
@@ -1680,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  */
@@ -1708,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)
@@ -1729,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;
 
@@ -1799,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++) {
@@ -2019,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];
 		}
@@ -2233,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);
@@ -2431,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);
@@ -2445,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;
@@ -2466,9 +2731,10 @@ MoleculeAddGaussianOrbitalShell(Molecule *mol, Int sym, Int nprims, Int a_idx)
 		case -1: shellp->sym = kGTOType_SP; shellp->ncomp = 4; break;
 		case 2:  shellp->sym = kGTOType_D;  shellp->ncomp = 6; break;
 		case -2: shellp->sym = kGTOType_D5; shellp->ncomp = 5; break;
-			/*  TODO: Support F/F7 type orbitals  */
-			/*	case 3: sp->sym = kGTOtype_F;  sp->ncomp = 10; break;
-			 case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break; */
+		case 3:  shellp->sym = kGTOType_F;  shellp->ncomp = 10; break;
+		case -3: shellp->sym = kGTOType_F7; shellp->ncomp = 7; break;
+		case 4:  shellp->sym = kGTOType_G;  shellp->ncomp = 15; break;
+		case -4: shellp->sym = kGTOType_G9; shellp->ncomp = 9; break;
 		default:
 			return -3;  /* Unsupported shell type  */
 	}
@@ -2481,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;
 }
 
@@ -2507,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)
 {
@@ -2536,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);
@@ -2549,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)
@@ -2565,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;
@@ -2708,6 +3090,7 @@ sSetupGaussianCoefficients(BasisSet *bset)
 					dp[3] = d * 1.425410941;
 					dp += 5;
 					break;
+				/*  TODO: Support F/F7 and G/G9 type orbitals  */
 			}
 		}
 	}
@@ -2728,7 +3111,7 @@ MoleculeLoadGaussianFchkFile(Molecule *mp, const char *fname, char **errbuf)
 	Int *iary;
 	Double *dary;
 	Atom *ap;
-	Vector *vp;
+/*	Vector *vp; */
 	Double w;
 
 	*errbuf = NULL;
@@ -2774,10 +3157,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) {
@@ -2860,13 +3244,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;
@@ -2926,9 +3307,10 @@ MoleculeLoadGaussianFchkFile(Molecule *mp, const char *fname, char **errbuf)
 					case -1: sp->sym = kGTOType_SP; sp->ncomp = 4; break;
 					case 2:  sp->sym = kGTOType_D;  sp->ncomp = 6; break;
 					case -2: sp->sym = kGTOType_D5; sp->ncomp = 5; break;
-						/*  TODO: Support F/F7 type orbitals  */
-						/*	case 3: sp->sym = kGTOtype_F;  sp->ncomp = 10; break;
-						 case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break; */
+					case 3:  sp->sym = kGTOType_F;  sp->ncomp = 10; break;
+					case -3: sp->sym = kGTOType_F7; sp->ncomp = 7; break;
+					case 4:  sp->sym = kGTOType_G;  sp->ncomp = 15; break;
+					case -4: sp->sym = kGTOType_G9; sp->ncomp = 9; break;
 					default:
 						s_append_asprintf(errbuf, "Line %d: unsupported shell type %d", lineNumber, iary[i]);
 						retval = 2;
@@ -3219,9 +3601,10 @@ MoleculeLoadGamessDatFile(Molecule *mol, const char *fname, char **errbuf)
 					vbuf[i].z = dval[3];
 				}
 				/*  Skip until a blank line is found  */
+				/*  2013.6.11. Line including "PM3" is also recognized as the end of atom  */
 				while ((status = sReadLineWithInterrupt(buf, sizeof buf, fp, &lineNumber)) > 0) {
 					for (j = 0; buf[j] == ' '; j++);
-					if (buf[j] == '\n')
+					if (buf[j] == '\n' || strncmp(buf + j, "PM3", 3) == 0)
 						break;
 				}
 				i++;
@@ -3282,7 +3665,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;
@@ -3317,7 +3700,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);
@@ -3362,7 +3745,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);
@@ -3884,8 +4267,9 @@ int
 MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 {
 	FILE *fp;
-	Int i, j, k, n1, n2, n3, n_aniso, nframes, nanchors;
+	Int i, j, k, n1, n2, n3, n_aniso, nframes, nanchors, n_uff;
 	Atom *ap;
+	char *p;
 	char bufs[6][8];
 
 	*errbuf = NULL;
@@ -3900,7 +4284,7 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 
 	fprintf(fp, "!:atoms\n");
 	fprintf(fp, "! idx seg_name res_seq res_name name type charge weight element atomic_number occupancy temp_factor int_charge\n");
-	n1 = n2 = n3 = n_aniso = nanchors = 0;
+	n1 = n2 = n3 = n_aniso = nanchors = n_uff = 0;
 	for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
 		strncpy(bufs[0], ap->segName, 4);
 		bufs[0][4] = 0;
@@ -3934,12 +4318,23 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 			n_aniso++;
 		if (ap->anchor != NULL)
 			nanchors++;
+		if (ap->uff_type[0] != 0)
+			n_uff++;
 		fprintf(fp, "%d %s %d %s %s %s %.5f %.5f %s %d %f %f %d\n", i, bufs[0], ap->resSeq, bufs[1], bufs[2], bufs[3], ap->charge, ap->weight, bufs[4], ap->atomicNumber, ap->occupancy, ap->tempFactor, ap->intCharge);
 	}
 	fprintf(fp, "\n");
 	
-	if (n1 > 0) {
-		fprintf(fp, "!:atoms_symop\n");
+	if (n_uff > 0) {
+		fprintf(fp, "!:uff_type\n");
+		fprintf(fp, "! idx uff_type\n");
+		for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
+			fprintf(fp, "%d %.5s\n", i, ap->uff_type);
+		}
+		fprintf(fp, "\n");
+	}
+	
+	if (n1 > 0) {
+		fprintf(fp, "!:atoms_symop\n");
 		fprintf(fp, "! idx symop symbase\n");
 		for (i = 0, ap = mp->atoms; i < mp->natoms; i++, ap = ATOM_NEXT(ap)) {
 			int n;
@@ -4019,6 +4414,15 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 		}
 		fprintf(fp, "\n");
 	}
+
+	if (mp->nbondOrders > 0) {
+		fprintf(fp, "!:bond_orders\n");
+		fprintf(fp, "! order1 order2 order3 order4\n");
+		for (i = 0; i < mp->nbondOrders; i++) {
+			fprintf(fp, "%.6f%c", mp->bondOrders[i], (i % 4 == 3 || i == mp->nbondOrders - 1 ? '\n' : ' '));
+		}
+		fprintf(fp, "\n");
+	}
 	
 	if (mp->nangles > 0) {
 		fprintf(fp, "!:angles\n");
@@ -4087,46 +4491,6 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 		fprintf(fp, "\n");
 	}
 	
-#if PIATOM
-	if (mp->npiatoms > 0) {
-		PiAtom *pp;
-		fprintf(fp, "!:pi_atoms\n");
-		fprintf(fp, "! name type n; a1 coeff1; a2 coeff2; ...; a_n coeff_n\n");
-		for (i = 0, pp = mp->piatoms; i < mp->npiatoms; i++, pp++) {
-			strncpy(bufs[0], pp->aname, 4);
-			bufs[0][4] = 0;
-			AtomTypeDecodeToString(pp->type, bufs[1]);
-			bufs[1][6] = 0;
-			for (j = 0; j < 2; j++) {
-				if (bufs[j][0] == 0) {
-					bufs[j][0] = '_';
-					bufs[j][1] = 0;
-				}
-			}
-			fprintf(fp, "%s %s %d\n", bufs[0], bufs[1], pp->connect.count);
-			ip = AtomConnectData(&pp->connect);
-			for (j = 0; j < pp->connect.count; j++) {
-				fprintf(fp, "%d %g\n", ip[j], (j < pp->ncoeffs ? pp->coeffs[j] : 0.0));
-			}
-		}
-		fprintf(fp, "\n");
-	}
-	
-	if (mp->npibonds > 0) {
-		fprintf(fp, "!:pi_atom_constructs\n");
-		fprintf(fp, "! a1 b1 c1 d1 a2 b2 c2 d2\n");
-		for (i = 0; i < mp->npibonds * 4; i++) {
-			j = mp->pibonds[i];
-			if (j >= ATOMS_MAX_NUMBER)
-				j = -2 - (j - ATOMS_MAX_NUMBER);
-			else if (j < 0)
-				j = -1;
-			fprintf(fp, "%d%c", j, (i % 8 == 7 || i == mp->npibonds * 4 - 1 ? '\n' : ' '));
-		}
-		fprintf(fp, "\n");
-	}
-#endif  /*  PIATOM  */
-
 	if (n_aniso > 0) {
 		fprintf(fp, "!:anisotropic_thermal_parameters\n");
 		fprintf(fp, "! b11 b22 b33 b12 b13 b23 [sigma; sb11 sb22 sb33 sb12 sb13 sb23]\n");
@@ -4162,6 +4526,7 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char **errbuf)
 		fprintf(fp, "cutoff %g\n", arena->cutoff);
 		fprintf(fp, "electro_cutoff %g\n", arena->electro_cutoff);
 		fprintf(fp, "pairlist_distance %g\n", arena->pairlist_distance);
+		fprintf(fp, "switch_distance %g\n", arena->switch_distance);
 		fprintf(fp, "temperature %g\n", arena->temperature);
 		fprintf(fp, "andersen_freq %d\n", arena->andersen_thermo_freq);
 		fprintf(fp, "andersen_coupling %g\n", arena->andersen_thermo_coupling);
@@ -4226,7 +4591,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");
@@ -4253,9 +4618,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;
 }
@@ -4364,24 +4861,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);
@@ -4833,58 +5312,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)
@@ -5092,6 +5519,14 @@ MoleculePrepareMDArena(Molecule *mol, int check_only, char **retmsg)
 		IntGroupRelease(ig3);
 	}
 	
+	{
+		/*  Update the path information of the molecule before MD setup  */
+		char *buf = (char *)malloc(4096);
+		MoleculeCallback_pathName(mol, buf, sizeof buf);
+		MoleculeSetPath(mol, buf);
+		free(buf);
+	}
+		
 	/*  Prepare parameters and internal information  */
 	msg = md_prepare(arena, check_only);
 	
@@ -5207,7 +5642,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);
@@ -5271,36 +5709,6 @@ MoleculeDeserialize(const char *data, Int length, Int *timep)
 				}
 				ptr += sizeof(Int) * (2 + n) + sizeof(Double) * n;
 			}
-#if PIATOM
-		} else if (strcmp(data, "PIATOM") == 0) {
-			const char *ptr2 = ptr + len;
-			mp->piatoms = NULL;
-			mp->npiatoms = 0;
-			while (ptr < ptr2) {
-				PiAtom pa;
-				memset(&pa, 0, sizeof(pa));
-				n = offsetof(PiAtom, connect);
-				memmove(&pa, ptr, n);
-				ptr += n;
-				n = *((Int *)ptr);
-				if (n > 0) {
-					AtomConnectResize(&pa.connect, n);
-					memmove(AtomConnectData(&pa.connect), ptr + sizeof(Int), n * sizeof(Int));
-				}
-				ptr += sizeof(Int) * (n + 1);
-				n = *((Int *)ptr);
-				if (n > 0) {
-					NewArray(&pa.coeffs, &pa.ncoeffs, sizeof(Double), n);
-					memmove(pa.coeffs, ptr + sizeof(Int), n * sizeof(Double));
-				}
-				ptr += sizeof(Int) + sizeof(Double) * n;
-				AssignArray(&mp->piatoms, &mp->npiatoms, sizeof(PiAtom), mp->npiatoms, &pa);
-			}
-		} else if (strcmp(data, "PIBOND") == 0) {
-			n = len / (sizeof(Int) * 4);
-			NewArray(&mp->pibonds, &mp->npibonds, sizeof(Int) * 4, n);
-			memmove(mp->pibonds, ptr, len);
-#endif  /*  PIATOM  */
 		} else if (strcmp(data, "TIME") == 0) {
 			if (timep != NULL)
 				*timep = *((Int *)ptr);
@@ -5586,58 +5994,6 @@ MoleculeSerialize(Molecule *mp, Int *outLength, Int *timep)
 		len_all += len;
 	}
 	
-#if PIATOM
-	/*  Pi-atoms  */
-	if (mp->npiatoms > 0) {
-		/*  Estimate the necessary storage first  */
-		len = 8 + sizeof(Int);
-		for (i = 0; i < mp->npiatoms; i++) {
-			len += offsetof(PiAtom, connect);  /*  Members before 'connect' is stored as they are  */
-			len += sizeof(Int) * (1 + mp->piatoms[i].connect.count); /* Array of Int's */
-			len += sizeof(Int) + sizeof(Double) * mp->piatoms[i].ncoeffs; /* Array of Double's */
-		}
-		ptr = (char *)realloc(ptr, len_all + len);
-		if (ptr == NULL)
-			goto out_of_memory;
-		p = ptr + len_all;
-		memmove(p, "PIATOM\0\0", 8);
-		*((Int *)(p + 8)) = len - (8 + sizeof(Int));
-		p += 8 + sizeof(Int);
-		for (i = 0; i < mp->npiatoms; i++) {
-			int len0;
-			PiAtom *pp = &(mp->piatoms[i]);
-			len0 = offsetof(PiAtom, connect);
-			memmove(p, pp, len0);
-			p += len0;
-			len0 = pp->connect.count * sizeof(Int);
-			*((Int *)p) = pp->connect.count;
-			if (len0 > 0)
-				memmove(p + sizeof(Int), AtomConnectData(&(pp->connect)), len0);
-			p += sizeof(Int) + len0;
-			len0 = pp->ncoeffs * sizeof(Double);
-			*((Int *)p) = pp->ncoeffs;
-			if (len0 > 0)
-				memmove(p + sizeof(Int), pp->coeffs, len0);
-			p += sizeof(Int) + len0;
-		}
-		len_all += len;
-	}
-	
-	/*  Pi-atom constructs  */
-	if (mp->npibonds > 0) {
-		len = 8 + sizeof(Int) + sizeof(Int) * 4 * mp->npibonds;
-		ptr = (char *)realloc(ptr, len_all + len);
-		if (ptr == NULL)
-			goto out_of_memory;
-		p = ptr + len_all;
-		memmove(p, "PIBOND\0\0", 8);
-		*((Int *)(p + 8)) = sizeof(Int) * 4 * mp->npibonds;
-		p += 8 + sizeof(Int);
-		memmove(p, mp->pibonds, sizeof(Int) * 4 * mp->npibonds);
-		len_all += len;
-	}
-#endif  /*  PIATOM  */
-
 	/*  Parameters  */
 	if (mp->par != NULL) {
 		int type;
@@ -5848,25 +6204,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)
@@ -5876,11 +6233,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;
@@ -5898,42 +6255,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;
@@ -6609,16 +6943,14 @@ MoleculeGetSymopForTransform(Molecule *mp, const Transform tf, Symop *symop, int
 	int i, j, n[3];
 	if (mp == NULL || mp->cell == NULL)
 		return -1;
-	if (mp->nsyms == 0)
-		return -2;
 	if (is_cartesian) {
 		TransformMul(t, tf, mp->cell->tr);
 		TransformMul(t, mp->cell->rtr, t);
 	} else {
 		memmove(t, tf, sizeof(Transform));
 	}
-	for (i = 0; i < mp->nsyms; i++) {
-		Transform *tp = mp->syms + i;
+	for (i = 0; i < mp->nsyms || i == 0; i++) {
+		Transform *tp = &(SYMMETRY_AT_INDEX(mp->syms, i));
 		for (j = 0; j < 9; j++) {
 			if (fabs((*tp)[j] - t[j]) > 1e-4)
 				break;
@@ -6672,24 +7004,28 @@ MoleculeTransformBySymop(Molecule *mp, const Vector *vpin, Vector *vpout, Symop
 	If indices is non-NULL, it should be an array of Int with at least 
 	IntGroupGetCount(group) entries, and on return it contains the
     indices of the expanded atoms (may be existing atoms if the expanded
-    atoms are already present)  */
+    atoms are already present)
+    If allowOverlap is non-zero, then the new atom is created even when the
+    coordinates coincide with the some other atom (special position) of the
+    same element; otherwise, such atom will not be created and the existing
+    atom is returned in indices[].  */
 int
-MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indices)
+MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indices, Int allowOverlap)
 {
-	int i, n, n0, n1, count, *table;
+	int i, n, n0, n1, n2, base, count, *table;
 	Atom *ap;
 	IntGroupIterator iter;
-	Transform tr;
-
+	Transform tr, t1;
+	Symop symop1;
+	Atom *ap2;
+	Vector nr, dr;
+	
 	if (mp == NULL || mp->natoms == 0 || group == NULL || (count = IntGroupGetCount(group)) == 0)
 		return -1;
-	if (symop.sym >= mp->nsyms)
+	if (symop.sym != 0 && symop.sym >= mp->nsyms)
 		return -2;
 
 	/*  Create atoms, with avoiding duplicates  */
-#if PIATOM
-	MoleculeInvalidatePiConnectionTable(mp);
-#endif
 	n0 = n1 = mp->natoms;
 	table = (int *)malloc(sizeof(int) * n0);
 	if (table == NULL)
@@ -6700,18 +7036,14 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice
 	MoleculeGetTransformForSymop(mp, symop, &tr, 0);
 	__MoleculeLock(mp);
 	for (i = 0; i < count; i++) {
-		int n2, base;
-		Symop symop1;
-		Atom *ap2;
-		Vector nr, dr;
 		n = IntGroupIteratorNext(&iter);
 		ap = ATOM_AT_INDEX(mp->atoms, n);
 		if (SYMOP_ALIVE(ap->symop)) {
 			/*  Calculate the cumulative symop  */
-			Transform t1;
+			Transform tr2;
 			MoleculeGetTransformForSymop(mp, ap->symop, &t1, 0);
-			TransformMul(t1, tr, t1);
-			if (MoleculeGetSymopForTransform(mp, t1, &symop1, 0) != 0) {
+			TransformMul(tr2, tr, t1);
+			if (MoleculeGetSymopForTransform(mp, tr2, &symop1, 0) != 0) {
 				if (indices != NULL)
 					indices[i] = -1;
 				continue;  /*  Skip this atom  */
@@ -6721,25 +7053,27 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice
 			symop1 = symop;
 			base = n;
 		}
+
+		/*  Calculate the expande position  */
+		MoleculeTransformBySymop(mp, &(ap->r), &nr, symop);
+		
 		/*  Is this expansion already present?  */
 		for (n2 = 0, ap2 = mp->atoms; n2 < n0; n2++, ap2 = ATOM_NEXT(ap2)) {
+			/*  Symmetry operation and the base atom are the same  */
 			if (ap2->symbase == base && SYMOP_EQUAL(symop1, ap2->symop))
 				break;
+			/*  Atomic number and the position are the same  */
+			if (ap2->atomicNumber == ap->atomicNumber && allowOverlap == 0) {
+				VecSub(dr, ap2->r, nr);
+				if (VecLength2(dr) < 1e-6)
+					break;
+			}
 		}
 		if (n2 < n0) {
 			/*  If yes, then skip it  */
 			if (indices != NULL)
 				indices[i] = n2;
 			continue;
-		}
-		/*  Is the expanded position coincides with itself?  */
-		MoleculeTransformBySymop(mp, &(ap->r), &nr, symop);
-		VecSub(dr, ap->r, nr);
-		if (VecLength2(dr) < 1e-6) {
-			/*  If yes, then this atom is included but no new atom is created  */
-			table[n] = n;
-			if (indices != NULL)
-				indices[i] = n;
 		} else {
 			/*  Create a new atom  */
 			Atom newAtom;
@@ -6748,9 +7082,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)
@@ -6761,21 +7095,35 @@ MoleculeAddExpandedAtoms(Molecule *mp, Symop symop, IntGroup *group, Int *indice
 	IntGroupIteratorRelease(&iter);
 
 	/*  Create bonds  */
-	for (i = 0; i < n0; i++) {
-		int b[2];
-		Int *cp;
-		b[0] = table[i];
-		if (b[0] < 0 || b[0] == i)
-			continue;
+	for (i = n0; i < n1; i++) {
+		Int b[2], j;
 		ap = ATOM_AT_INDEX(mp->atoms, i);
-		cp = AtomConnectData(&ap->connect);
-		for (n = 0; n < ap->connect.count; n++) {
-			b[1] = table[cp[n]];
-			if (b[1] < 0)
-				continue;
-			if (b[1] > n0 && b[0] > b[1])
-				continue;
-			MoleculeAddBonds(mp, 1, b, NULL, 1);
+		if (SYMOP_ALIVE(ap->symop) && MoleculeGetTransformForSymop(mp, ap->symop, &tr, 1) == 0) {
+			/*  For each connected atom, look for the transformed atom  */
+			Int *cp;
+			ap2 = ATOM_AT_INDEX(mp->atoms, ap->symbase);
+			cp = AtomConnectData(&ap2->connect);
+			n2 = ap2->connect.count;
+			for (n = 0; n < n2; n++) {
+				Atom *apn = ATOM_AT_INDEX(mp->atoms, cp[n]);
+				nr = apn->r;
+				TransformVec(&nr, tr, &nr);
+				/*  Look for the bonded atom transformed by ap->symop  */
+				for (j = 0, ap2 = mp->atoms; j < mp->natoms; j++, ap2 = ATOM_NEXT(ap2)) {
+					if (ap2->symbase == cp[n] && SYMOP_EQUAL(ap->symop, ap2->symop))
+						break;
+					VecSub(dr, nr, ap2->r);
+					if (ap2->atomicNumber == apn->atomicNumber && VecLength2(dr) < 1e-6)
+						break;
+				}
+				if (j < mp->natoms) {
+					/*  Bond i-j is created  */
+					b[0] = i;
+					b[1] = j;
+					if (MoleculeLookupBond(mp, b[0], b[1]) < 0)
+						MoleculeAddBonds(mp, 1, b, NULL, 1);
+				}
+			}
 		}
 	}
 	mp->needsMDRebuild = 1;
@@ -6875,8 +7223,10 @@ MoleculeAmendBySymmetry(Molecule *mp, IntGroup *group, IntGroup **groupout, Vect
 			free(vp);
 		}
 	} else {
-		*groupout = NULL;
-		*vpout = NULL;
+		if (groupout != NULL && vpout != NULL) {
+			*groupout = NULL;
+			*vpout = NULL;
+		}
 	}
 	return count;
 }
@@ -6998,7 +7348,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  */
@@ -7079,11 +7431,18 @@ MoleculeCreateAnAtom(Molecule *mp, const Atom *ap, int pos)
 		for (i = 0, api = ATOM_AT_INDEX(mp->atoms, i); i < mp->natoms; i++, api = ATOM_NEXT(api)) {
 			int j;
 			Int *cp;
+			cp = AtomConnectData(&api->connect);
 			for (j = 0; j < api->connect.count; j++) {
-				cp = AtomConnectData(&api->connect);
 				if (cp[j] >= pos)
 					cp[j]++;
 			}
+			if (api->anchor != NULL) {
+				cp = AtomConnectData(&api->anchor->connect);
+				for (j = 0; j < api->anchor->connect.count; j++) {
+					if (cp[j] >= pos)
+						cp[j]++;
+				}
+			}
 		}
 		for (i = 0; i < mp->nbonds * 2; i++) {
 			if (mp->bonds[i] >= pos)
@@ -7101,9 +7460,6 @@ MoleculeCreateAnAtom(Molecule *mp, const Atom *ap, int pos)
 			if (mp->impropers[i] >= pos)
 				mp->impropers[i]++;
 		}
-#if PIATOM
-		MoleculeInvalidatePiConnectionTable(mp);
-#endif
 	}
 	mp->nframes = -1;  /*  Should be recalculated later  */
 	MoleculeIncrementModifyCount(mp);
@@ -7229,9 +7585,7 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
 	act = NULL;
 
 	__MoleculeLock(dst);
-#if PIATOM
-	MoleculeInvalidatePiConnectionTable(dst);
-#endif
+
 	nsrc = src->natoms;
 	ndst = dst->natoms;
 	if (resSeqOffset < 0)
@@ -7357,6 +7711,17 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
 				goto panic;
 			/*  Copy the items  */
 			memmove(*items + n1 * nsize, *items_src, sizeof(Int) * nsize * n2);
+			if (i == 0) {
+				/*  Copy the bond order info if present */
+				Int nn1 = dst->nbondOrders;
+				if (dst->bondOrders != NULL || src->bondOrders != NULL) {
+					if (AssignArray(&dst->bondOrders, &dst->nbondOrders, sizeof(Double), dst->nbonds - 1, NULL) == NULL)
+						goto panic;
+					memset(dst->bondOrders + nn1, 0, sizeof(Double) * (dst->nbonds - nn1));
+					if (src->bondOrders != NULL)
+						memmove(dst->bondOrders + n1, src->bondOrders, sizeof(Double) * n2);
+				}
+			}
 		}
 		/*  Renumber  */
 		for (j = 0; j < n1 * nsize; j++)
@@ -7487,81 +7852,6 @@ MoleculeMerge(Molecule *dst, Molecule *src, IntGroup *where, Int resSeqOffset, I
 		}
 	}
 
-#if PIATOM
-	/*  Renumber the existing pi-atoms  */
-	if (dst->npiatoms > 0) {
-		for (i = 0; i < dst->npiatoms; i++) {
-			PiAtom *pp;
-			pp = &dst->piatoms[i];
-			cp = AtomConnectData(&pp->connect);
-			for (j = 0; j < pp->connect.count; j++) {
-				cp[j] = old2new[cp[j]];
-			}
-		}
-		if (dst->npibonds > 0) {
-			cp = dst->pibonds;
-			for (i = 0; i < dst->npibonds * 4; i++) {
-				if (cp[i] < 0)
-					continue;
-				else if (cp[i] >= ATOMS_MAX_NUMBER)
-					continue;
-				else {
-					cp[i] = old2new[cp[i]];
-				}
-			}
-		}
-	}
-	
-	/*  Copy the pi-atoms  */
-	if (src->npiatoms > 0 && forUndo == 0) {
-		int nsrcp, ndstp;
-		nsrcp = src->npiatoms;
-		ndstp = dst->npiatoms;
-		if (AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), nsrcp + ndstp - 1, NULL) == NULL)
-			goto panic;
-		for (i = 0; i < nsrcp; i++) {
-			PiAtom *pp;
-			pp = &dst->piatoms[ndstp + i];
-			PiAtomDuplicate(pp, &src->piatoms[i]);
-			cp = AtomConnectData(&pp->connect);
-			for (j = 0; j < pp->connect.count; j++) {
-				cp[j] = old2new[ndst + cp[j]];
-			}
-			if (nactions != NULL) {
-				/*  This is very inefficient, yet should not cause big problems
-				    because the number of piatoms in the molecule is usually limited.  */
-				act = MolActionNew(gMolActionRemoveOnePiAtom, ndstp + i);
-				AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-				act = NULL;
-			}
-		}
-		if (src->npibonds > 0) {
-			n1 = src->npibonds;
-			n2 = dst->npibonds;
-			if (AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, n1 + n2 - 1, NULL) == NULL)
-				goto panic;
-			cp = &dst->pibonds[n2 * 4];
-			memmove(cp, src->pibonds, sizeof(Int) * 4 * n1);
-			for (i = 0; i < 4 * n1; i++) {
-				/*  Renumber the pi-atom constructs  */
-				if (cp[i] < 0)
-					continue;
-				else if (cp[i] < ATOMS_MAX_NUMBER)
-					cp[i] = old2new[ndst + cp[i]];  /*  Ordinary atoms  */
-				else
-					cp[i] += ndstp;  /*  pi-atoms  */
-			}
-			if (nactions != NULL) {
-				ig = IntGroupNewWithPoints(n2, n1, -1);
-				act = MolActionNew(gMolActionRemovePiBonds, ig);
-				AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-				act = NULL;
-				IntGroupRelease(ig);
-			}
-		}
-	}
-#endif  /*  PIATOM  */
-
 	MoleculeCleanUpResidueTable(dst);
 	
 	free(new2old);
@@ -7606,9 +7896,6 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
 		return 1;  /*  Bad parameter  */
 
 	__MoleculeLock(src);
-#if PIATOM
-	MoleculeInvalidatePiConnectionTable(src);
-#endif
 	
 	act = NULL;
 	
@@ -7881,14 +8168,26 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
 					goto panic;
 				if (sCopyElementsFromArrayAtPositions(*items, *nitems, *items_dst, sizeof(Int) * nsize, move_g) != 0)
 					goto panic;
+				if (i == 0 && src->bondOrders != NULL) {
+					if (AssignArray(&dst->bondOrders, &dst->nbondOrders, sizeof(Double), n3 - 1, NULL) == NULL)
+						goto panic;
+					if (sCopyElementsFromArrayAtPositions(src->bondOrders, src->nbondOrders, dst->bondOrders, sizeof(Double), move_g) != 0)
+						goto panic;
+				}
 			}
 			/*  Remove from src  */
 			if (moveFlag && forUndo == 0) {
 				if (nactions != NULL) {
 					Int k, *ip;
+					Double *dp;
 					ip = (Int *)malloc(sizeof(Int) * nsize * n2);
 					for (j = 0; (k = IntGroupGetNthPoint(del_g, j)) >= 0; j++)
 						memmove(ip + j * nsize, *items + k * nsize, sizeof(Int) * nsize);
+					if (i == 0 && src->bondOrders != NULL) {
+						dp = (Double *)malloc(sizeof(Double) * n2);
+						for (j = 0; (k = IntGroupGetNthPoint(del_g, j)) >= 0; j++)
+							dp[j] = src->bondOrders[k];
+					} else dp = NULL;
 					switch (i) {
 						case 0:
 							act = MolActionNew(gMolActionAddBondsForUndo, n2 * nsize, ip, del_g); break;
@@ -7904,6 +8203,12 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
 						act = NULL;
 					}
 					free(ip);
+					if (dp != NULL) {
+						act = MolActionNew(gMolActionAssignBondOrders, n2, dp, del_g);
+						AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
+						act = NULL;
+						free(dp);
+					}
 				}
 				if (sRemoveElementsFromArrayAtPositions(*items, *nitems, NULL, sizeof(Int) * nsize, del_g) != 0)
 					goto panic;
@@ -8000,149 +8305,6 @@ sMoleculeUnmergeSub(Molecule *src, Molecule **dstp, IntGroup *where, int resSeqO
 		}
 	}
 
-#if PIATOM
-	/*  Copy the pi-atoms  */
-	if (src->npiatoms > 0) {
-		PiAtom *pp1, *pp2;
-		Int *patoms_old2new;
-		n1 = src->npiatoms;
-		patoms_old2new = (Int *)calloc(sizeof(Int), src->npiatoms);
-		if (patoms_old2new == NULL)
-			goto panic;
-		for (i = 0, pp1 = src->piatoms; i < src->npiatoms; i++, pp1++) {
-			/*  Is this entry to be copied to dst?  */
-			cp = AtomConnectData(&pp1->connect);
-			for (j = pp1->connect.count - 1; j >= 0; j--) {
-				if (old2new[cp[j]] < nsrcnew)
-					break;
-			}
-			if (j < 0) {
-				/*  Copy this entry  */
-				patoms_old2new[i] = dst->npiatoms;
-				pp2 = AssignArray(&dst->piatoms, &dst->npiatoms, sizeof(PiAtom), dst->npiatoms, NULL);
-				PiAtomDuplicate(pp2, pp1);
-				cp = AtomConnectData(&pp2->connect);
-				for (j = 0; j < pp2->connect.count; j++)
-					cp[j] = old2new[cp[j]] - nsrcnew;
-			} else {
-				patoms_old2new[i] = -1;
-			}
-		}
-		/*  Copy the piatom constructs to dst  */
-		for (i = 0; i < src->npibonds; i++) {
-			for (j = 0; j < 4; j++) {
-				n2 = src->pibonds[i * 4 + j];
-				if (n2 >= ATOMS_MAX_NUMBER) {
-					if (patoms_old2new[n2 - ATOMS_MAX_NUMBER] < 0)
-						break;
-				} else if (n2 >= 0) {
-					if (old2new[n2] < nsrcnew)
-						break;
-				}
-			}
-			if (j >= 4) {
-				/*  Copy this entry  */
-				cp = (Int *)AssignArray(&dst->pibonds, &dst->npibonds, sizeof(Int) * 4, dst->npibonds, &src->pibonds[i * 4]);
-				for (j = 0; j < 4; j++) {
-					if (cp[j] >= ATOMS_MAX_NUMBER) {
-						cp[j] = patoms_old2new[cp[j] - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
-					} else if (cp[j] >= 0) {
-						cp[j] = old2new[cp[j]] - nsrcnew;
-					}
-				}
-			}
-		}
-		if (moveFlag) {
-			Int npibonds_to_go, *pibonds_to_go;
-			IntGroup *pibonds_group;
-
-			/*  Remove the piatom entries containing non-remaining atoms. Note: the piatom
-			 entries that do not remain in src and not copied to dst will disappear.  */
-			n2 = 0;
-			for (i = 0; i < src->npiatoms; i++) {
-				/*  Is this entry to be removed?  */
-				pp1 = &src->piatoms[i];
-				cp = AtomConnectData(&pp1->connect);
-				for (j = pp1->connect.count - 1; j >= 0; j--) {
-					if (old2new[cp[j]] >= nsrcnew)
-						break;
-				}
-				patoms_old2new[i] = (j < 0 ? n2++ : -1);
-			}
-
-			/*  Remove pibonds first (necessary for undo)  */
-			npibonds_to_go = 0;
-			pibonds_to_go = NULL;
-			pibonds_group = NULL;
-			for (i = src->npibonds - 1; i >= 0; i--) {
-				for (j = 0; j < 4; j++) {
-					n3 = src->pibonds[i * 4 + j];
-					if (n3 >= ATOMS_MAX_NUMBER) {
-						if (patoms_old2new[n3 - ATOMS_MAX_NUMBER] < 0)
-							break;
-					} else if (n3 >= 0) {
-						if (old2new[n3] >= nsrcnew)
-							break;
-					}
-				}
-				if (j < 4) {
-					/*  Remove  */
-					if (nactions != NULL) {
-						/*  Since we are scanning pibonds[] from the end, the new entry should be inserted
-						    at the top of the pibonds_to_go[] array.  */
-						InsertArray(&pibonds_to_go, &npibonds_to_go, sizeof(Int) * 4, 0, 1, src->pibonds + i * 4);
-						if (pibonds_group == NULL)
-							pibonds_group = IntGroupNew();
-						IntGroupAdd(pibonds_group, i, 1);
-					}
-					DeleteArray(&src->pibonds, &src->npibonds, sizeof(Int) * 4, i, 1, NULL);
-				} else {
-					/*  Renumber  */
-					cp = &src->pibonds[i * 4];
-					for (j = 0; j < 4; j++) {
-						n3 = cp[j];
-						if (n3 >= ATOMS_MAX_NUMBER)
-							cp[j] = patoms_old2new[n3 - ATOMS_MAX_NUMBER] + ATOMS_MAX_NUMBER;
-						else if (n3 >= 0)
-							cp[j] = old2new[n3];
-					}
-				}
-			}
-			if (nactions != NULL && pibonds_to_go != NULL) {
-				act = MolActionNew(gMolActionInsertPiBonds, pibonds_group, npibonds_to_go * 4, pibonds_to_go);
-				AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-				act = NULL;
-				free(pibonds_to_go);
-			}
-			if (pibonds_group != NULL)
-				IntGroupRelease(pibonds_group);
-
-			for (i = src->npiatoms - 1; i >= 0; i--) {
-				pp1 = &src->piatoms[i];
-				if (patoms_old2new[i] < 0) {
-					/*  Remove the entries  */
-					/*  (If forUndo is true, these entries should already have been removed.)  */
-					if (nactions != NULL) {
-						act = MolActionNew(gMolActionInsertOnePiAtom, i, 4, pp1->aname, pp1->type, pp1->connect.count, AtomConnectData(&pp1->connect), pp1->ncoeffs, pp1->coeffs);
-						AssignArray(actions, nactions, sizeof(MolAction *), *nactions, &act);
-						act = NULL;
-					}
-					PiAtomClean(pp1);
-					DeleteArray(&src->piatoms, &src->npiatoms, sizeof(PiAtom), i, 1, NULL);
-				} else {
-					/*  Renumber  */
-					cp = AtomConnectData(&pp1->connect);
-					for (j = 0; j < pp1->connect.count; j++) {
-						cp[j] = old2new[cp[j]];
-					}
-				}
-			}
-		}
-		
-		free(patoms_old2new);
-	}
-#endif  /*  PIATOM  */
-	
 	/*  Clean up  */
 	IntGroupRelease(remain_g);
 	MoleculeCleanUpResidueTable(src);
@@ -8247,10 +8409,13 @@ MoleculeExtract(Molecule *src, Molecule **dstp, IntGroup *where, int dummyFlag)
 int
 MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, Int autoGenerate)
 {
-	int i, j, n1, n2;
-	Atom *ap;
-	Int *cp;
-
+	Int nangles, ndihedrals;
+	Int *angles, *dihedrals;
+	Int i, j, k, kk, n1, n2, cn1, cn2;
+	Int *cp1, *cp2;
+	Int temp[4];
+	Atom *ap1, *ap2, *ap3;
+	
 	if (mp == NULL || bonds == NULL || nbonds <= 0)
 		return 0;
 	if (mp->noModifyTopology)
@@ -8266,37 +8431,39 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, In
 		__MoleculeUnlock(mp);
 		return -4;  /*  Out of memory  */
 	}
+	if (mp->bondOrders != NULL) {
+		/*  Expand the bond order info (all new entries are zero)  */
+		Double *dp = (Double *)calloc(sizeof(Double), nbonds);
+		if (dp == NULL)
+			return -4;
+		if (AssignArray(&(mp->bondOrders), &(mp->nbondOrders), sizeof(Double), n1 + nbonds - 1, NULL) == NULL
+			|| sInsertElementsToArrayAtPositions(mp->bondOrders, n1, dp, nbonds, sizeof(Double), where) != 0) {
+			__MoleculeUnlock(mp);
+			free(dp);
+			return -4;
+		}
+		free(dp);
+	}
 	
-	/*  Add connects[]  */
+	angles = dihedrals = NULL;
+	nangles = ndihedrals = 0;
+	
+	/*  Add connects[], and angles/dihedrals (if autoGenerate is true)  */
 	for (i = 0; i < nbonds; i++) {
+		
+		/*  One entry at time  */
+		/*  (Otherwise, duplicate entries of angles and dihedrals result)  */
 		n1 = bonds[i * 2];
 		n2 = bonds[i * 2 + 1];
-		ap = ATOM_AT_INDEX(mp->atoms, n1);
-		cp = AtomConnectData(&ap->connect);
-		AtomConnectInsertEntry(&ap->connect, -1, n2);
-		ap = ATOM_AT_INDEX(mp->atoms, n2);
-		cp = AtomConnectData(&ap->connect);
-		AtomConnectInsertEntry(&ap->connect, -1, n1);
-	}
+		
+		ap1 = ATOM_AT_INDEX(mp->atoms, n1);
+		AtomConnectInsertEntry(&ap1->connect, -1, n2);
+		ap2 = ATOM_AT_INDEX(mp->atoms, n2);
+		AtomConnectInsertEntry(&ap2->connect, -1, n1);
 	
-	/*  Add angles, dihedrals, impropers  */
-	if (autoGenerate) {
-		Int nangles, ndihedrals;
-		Int *angles, *dihedrals;
-		Int k, kk;
-		Int *cp1, *cp2;
-		Int temp[4];
-		Atom *ap1, *ap2, *ap3;
-
-		angles = dihedrals = NULL;
-		nangles = ndihedrals = 0;
-
-		for (i = 0; i < nbonds; i++) {
+		/*  Add angles and dihedrals  */
+		if (autoGenerate) {
 			AtomConnect *ac1, *ac2;
-			n1 = bonds[i * 2];
-			n2 = bonds[i * 2 + 1];
-			ap1 = ATOM_AT_INDEX(mp->atoms, n1);
-			ap2 = ATOM_AT_INDEX(mp->atoms, n2);
 			if (ap1->anchor == NULL || ap2->anchor == NULL) {
 				/*  N1-N2-{XY} or N2-N1-{XY} angles (X: connected atom, Y: constitute atom of pi-anchor)  */
 				for (j = 0; j < 4; j++) {
@@ -8307,16 +8474,22 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, In
 						case 3: if (ap1->anchor == NULL) continue; else ac1 = &ap1->anchor->connect; break; /* N2-N1-Y */
 					}
 					cp1 = AtomConnectData(ac1);
-					for (k = 0; k < ac1->count; k++) {
+					cn1 = ac1->count;
+					for (k = 0; k < cn1; k++) {
 						temp[2] = cp1[k];
 						if (temp[2] == temp[0])
 							continue;
+						ap3 = ATOM_AT_INDEX(mp->atoms, temp[2]);
+						if (ap3->anchor != NULL) {
+							/*  Avoid X-anchor-anchor angle (anchor-X-anchor is allowed)  */
+							if ((j < 2 && ap2->anchor != NULL) || (j >= 2 && ap1->anchor != NULL))
+								continue;
+						}
 						if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
 							goto panic;
-						/*  Dihedrals N1-N2-X-X or N2-N1-X-X  */
+						/*  Dihedrals N1-N2-X-{XY} or N2-N1-X-{XY}  */
 						if (j == 1 || j == 3)
 							continue;
-						ap3 = ATOM_AT_INDEX(mp->atoms, temp[2]);
 						cp2 = AtomConnectData(&ap3->connect);
 						for (kk = 0; kk < ap3->connect.count; kk++) {
 							temp[3] = cp2[kk];
@@ -8325,52 +8498,70 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, In
 							if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
 								goto panic;
 						}
+						if (ap3->anchor != NULL) {
+							/*  N1-N2-X-Y or N2-N1-X-Y  */
+							/*  for Y, only the first constitute atom is considered  */
+							cp2 = AtomConnectData(&ap3->anchor->connect);
+							temp[3] = cp2[0];
+							if (temp[3] == temp[0] || temp[3] == temp[1])
+								continue;
+							if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+								goto panic;
+						}
 					}
 				}
 			}
-			/*  X-N1-N2-X angles  */
-			if (ap1->anchor == NULL)
+			/*  X-N1-N2-X dihedrals  */
+			/*  Y-N1-N2-anchor is allowed, but the force may be zero if the angle N1-N2-anchor is */
+			/*  close to 180 deg (e.g. in ferrocene, C-anchor-Fe-anchor dihedral should be k=0)  */
+			if (ap1->anchor == NULL) {
 				ac1 = &ap1->connect;
-			else ac1 = &ap1->anchor->connect;
-			if (ap2->anchor == NULL)
+				cn1 = ac1->count;
+			} else {
+				ac1 = &ap1->anchor->connect;
+				cn1 = 1;  /*  Only the first constitute atom of pi-anchor is considered  */
+			}
+			if (ap2->anchor == NULL) {
 				ac2 = &ap2->connect;
-			else ac2 = &ap2->anchor->connect;
+				cn2 = ac2->count;
+			} else {
+				ac2 = &ap2->anchor->connect;
+				cn2 = 1;  /*  Only the first constitute atom of pi-anchor is considered  */
+			}
 			temp[1] = n1;
 			temp[2] = n2;
 			cp1 = AtomConnectData(ac1);
 			cp2 = AtomConnectData(ac2);
-			for (j = 0; j < ac1->count; j++) {
+			for (j = 0; j < cn1; j++) {
 				temp[0] = cp1[j];
 				if (temp[0] == temp[2])
 					continue;
-				if (ATOM_AT_INDEX(mp->atoms, temp[0])->anchor != NULL)
-					continue;
-				for (k = 0; k < ac2->count; k++) {
+				for (k = 0; k < cn2; k++) {
 					temp[3] = cp2[k];
 					if (temp[3] == temp[0] || temp[3] == temp[1])
 						continue;
-					if (ATOM_AT_INDEX(mp->atoms, temp[3])->anchor != NULL)
-						continue;
 					if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
 						goto panic;
 				}
 			}
 		}
-		temp[0] = kInvalidIndex;  /*  For termination  */
-		if (angles != NULL) {
-			if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
-				goto panic;
-			MoleculeAddAngles(mp, angles, NULL);
-			free(angles);
-		}
-		if (dihedrals != NULL) {
-			if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
-				goto panic;
-			MoleculeAddDihedrals(mp, dihedrals, NULL);
-			free(dihedrals);
-		}
 	}
 	
+	if (angles != NULL) {
+		temp[0] = kInvalidIndex;
+		if (AssignArray(&angles, &nangles, sizeof(Int) * 3, nangles, temp) == NULL)
+			goto panic;
+		MoleculeAddAngles(mp, angles, NULL);
+		free(angles);
+	}
+	if (dihedrals != NULL) {
+		temp[0] = kInvalidIndex;
+		if (AssignArray(&dihedrals, &ndihedrals, sizeof(Int) * 4, ndihedrals, temp) == NULL)
+			goto panic;
+		MoleculeAddDihedrals(mp, dihedrals, NULL);
+		free(dihedrals);
+	}
+
 	MoleculeIncrementModifyCount(mp);
 	mp->needsMDRebuild = 1;
 	__MoleculeUnlock(mp);
@@ -8397,8 +8588,8 @@ MoleculeAddBonds(Molecule *mp, Int nbonds, const Int *bonds, IntGroup *where, In
 int
 MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved, IntGroup **outRemovedPos)
 {
-	Int i, j, n1, n2;
-	Int *ip, na, nd, ni;
+	Int i, j, n1, n2, nw;
+	Int *ip, *jp, na, nd, ni;
 	IntGroup *ag, *dg, *ig;
 	Atom *ap;
 	IntGroupIterator iter;
@@ -8432,29 +8623,41 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 			}
 		}
 	}
-	IntGroupIteratorReset(&iter);
 	
 	/*  Remove bonds, angles, dihedrals, impropers  */
-	ag = dg = ig = NULL;
+	ag = IntGroupNew();
+	dg = ig = NULL;
 	na = nd = ni = 0;
+	
+	nw = IntGroupGetCount(where);
+	jp = (Int *)malloc(sizeof(Int) * nw * 2);
+	j = 0;
+	IntGroupIteratorReset(&iter);
 	while ((i = IntGroupIteratorNext(&iter)) >= 0) {
-		n1 = mp->bonds[i * 2];
-		n2 = mp->bonds[i * 2 + 1];
-		for (j = 0; j < mp->nangles; j++) {
-			ip = mp->angles + j * 3;
+		jp[j++] = mp->bonds[i * 2];
+		jp[j++] = mp->bonds[i * 2 + 1];
+	}
+	IntGroupIteratorRelease(&iter);
+
+	for (i = 0, ip = mp->angles; i < mp->nangles; i++, ip += 3) {
+		for (j = 0; j < nw; j++) {
+			n1 = jp[j * 2];
+			n2 = jp[j * 2 + 1];
 			if ((ip[0] == n1 && ip[1] == n2)
-			 || (ip[1] == n1 && ip[0] == n2)
-			 || (ip[1] == n1 && ip[2] == n2)
-			 || (ip[2] == n1 && ip[1] == n2)) {
-				if (ag == NULL)
-					ag = IntGroupNew();
-				if (IntGroupAdd(ag, j, 1) != 0)
+				|| (ip[1] == n1 && ip[0] == n2)
+				|| (ip[1] == n1 && ip[2] == n2)
+				|| (ip[2] == n1 && ip[1] == n2)) {
+				if (IntGroupAdd(ag, i, 1) != 0)
 					goto panic;
 				na++;
+				break;
 			}
 		}
-		for (j = 0; j < mp->ndihedrals; j++) {
-			ip = mp->dihedrals + j * 4;
+	}
+	for (i = 0, ip = mp->dihedrals; i < mp->ndihedrals; i++, ip += 4) {
+		for (j = 0; j < nw; j++) {
+			n1 = jp[j * 2];
+			n2 = jp[j * 2 + 1];
 			if ((ip[0] == n1 && ip[1] == n2)
 			 || (ip[1] == n1 && ip[0] == n2)
 			 || (ip[1] == n1 && ip[2] == n2)
@@ -8463,13 +8666,17 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 			 || (ip[3] == n1 && ip[2] == n2)) {
 				if (dg == NULL)
 					dg = IntGroupNew();
-				if (IntGroupAdd(dg, j, 1) != 0)
+				if (IntGroupAdd(dg, i, 1) != 0)
 					goto panic;
 				nd++;
+				break;
 			}
 		}
-		for (j = 0; j < mp->nimpropers; j++) {
-			ip = mp->impropers + j * 4;
+	}
+	for (i = 0, ip = mp->impropers; i < mp->nimpropers; i++, ip += 4) {
+		for (j = 0; j < nw; j++) {
+			n1 = jp[j * 2];
+			n2 = jp[j * 2 + 1];
 			if ((ip[0] == n1 && ip[2] == n2)
 			 || (ip[1] == n1 && ip[2] == n2)
 			 || (ip[3] == n1 && ip[2] == n2)
@@ -8478,14 +8685,15 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 			 || (ip[3] == n2 && ip[2] == n1)) {
 				if (ig == NULL)
 					ig = IntGroupNew();
-				if (IntGroupAdd(ig, j, 1) != 0)
+				if (IntGroupAdd(ig, i, 1) != 0)
 					goto panic;
 				ni++;
+				break;
 			}
 		}
 	}
-	IntGroupIteratorRelease(&iter);
-
+	free(jp);
+	
 	if (sRemoveElementsFromArrayAtPositions(mp->bonds, mp->nbonds, bonds, sizeof(Int) * 2, where) != 0)
 		goto panic;
 	mp->nbonds -= IntGroupGetCount(where);
@@ -8493,6 +8701,15 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 		free(mp->bonds);
 		mp->bonds = NULL;
 	}
+	if (mp->bondOrders != NULL) {
+		if (sRemoveElementsFromArrayAtPositions(mp->bondOrders, mp->nbondOrders, NULL, sizeof(Double), where) != 0)
+			goto panic;
+		mp->nbondOrders -= IntGroupGetCount(where);
+		if (mp->nbondOrders == 0) {
+			free(mp->bondOrders);
+			mp->bondOrders = NULL;
+		}
+	}
 	if (na == 0 && nd == 0 && ni == 0)
 		ip = NULL;
 	else
@@ -8512,6 +8729,11 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 		IntGroupRelease(ig);
 	}
 
+	if (IntGroupGetCount(ag) == 0) {
+		IntGroupRelease(ag);
+		ag = NULL;
+	}
+	
 	*outRemoved = ip;
 	*outRemovedPos = ag;
 
@@ -8528,6 +8750,55 @@ MoleculeDeleteBonds(Molecule *mp, Int *bonds, IntGroup *where, Int **outRemoved,
 }
 
 int
+MoleculeAssignBondOrders(Molecule *mp, const Double *orders, IntGroup *where)
+{
+	Int i, j;
+	IntGroupIterator iter;
+	if (mp == NULL || orders == NULL || mp->nbonds == 0)
+		return 0;
+	if (mp->noModifyTopology)
+		return -4;  /*  Prohibited operation  */
+	if (mp->bondOrders == NULL) {
+		AssignArray(&mp->bondOrders, &mp->nbondOrders, sizeof(Double), mp->nbonds - 1, NULL);
+		memset(mp->bondOrders, 0, sizeof(Double) * mp->nbondOrders);
+	}
+	IntGroupIteratorInit(where, &iter);
+	j = 0;
+	while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+		if (i >= mp->nbondOrders)
+			break;
+		mp->bondOrders[i] = orders[j++];
+	}
+	IntGroupIteratorRelease(&iter);
+	return 0;
+}
+
+int
+MoleculeGetBondOrders(Molecule *mp, Double *outOrders, IntGroup *where)
+{
+	Int i, j;
+	IntGroupIterator iter;
+	if (mp == NULL || mp->nbonds == 0)
+		return 0;
+	if (mp->bondOrders == NULL) {
+		/*  Returns all zero  */
+		i = IntGroupGetCount(where);
+		for (j = 0; j < i; j++)
+			outOrders[j] = 0.0;
+	} else {
+		IntGroupIteratorInit(where, &iter);
+		j = 0;
+		while ((i = IntGroupIteratorNext(&iter)) >= 0) {
+			if (i < mp->nbondOrders)
+				outOrders[j] = mp->bondOrders[i];
+			else outOrders[j] = 0.0;
+			j++;
+		}
+	}
+	return 0;
+}
+
+int
 MoleculeAddAngles(Molecule *mp, const Int *angles, IntGroup *where)
 {
 	int n1, nc;
@@ -8568,7 +8839,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) {
@@ -8620,7 +8891,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) {
@@ -8672,7 +8943,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) {
@@ -9117,12 +9388,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);
@@ -9163,7 +9434,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,
@@ -9392,10 +9663,6 @@ sMoleculeReorder(Molecule *mp)
 	free(newAtoms);
 	free(old2new);
 	free(apArray);
-	
-#if PIATOM
-	MoleculeInvalidatePiConnectionTable(mp);
-#endif
 }
 
 /*  Renumber atoms  */
@@ -9456,12 +9723,19 @@ MoleculeRenumberAtoms(Molecule *mp, const Int *new2old, Int *old2new_out, Int is
 	for (i = 0; i < mp->nimpropers * 4; i++) {
 		mp->impropers[i] = old2new[mp->impropers[i]];
 	}
+	/*  Renumber the connection table and pi anchor table  */
 	for (i = 0; i < mp->natoms; i++) {
 		Atom *ap = ATOM_AT_INDEX(saveAtoms, i);
 		Int *ip = AtomConnectData(&ap->connect);
 		for (j = 0; j < ap->connect.count; j++, ip++)
 			*ip = old2new[*ip];
+		if (ap->anchor != NULL) {
+			ip = AtomConnectData(&ap->anchor->connect);
+			for (j = 0; j < ap->anchor->connect.count; j++, ip++)
+				*ip = old2new[*ip];
+		}
 	}
+	
 	if (mp->par != NULL) {
 		/*  Renumber the parameters  */
 		int n;
@@ -9475,35 +9749,11 @@ MoleculeRenumberAtoms(Molecule *mp, const Int *new2old, Int *old2new_out, Int is
 		}
 	}
 	
-#if PIATOM
-	if (mp->npiatoms) {
-		/*  Renumber the pi-atoms  */
-		for (i = 0; i < mp->npiatoms; i++) {
-			PiAtom *pp = &mp->piatoms[i];
-			Int *cp = AtomConnectData(&pp->connect);
-			for (j = 0; j < pp->connect.count; j++)
-				cp[j] = old2new[cp[j]];
-		}
-	}
-	
-	if (mp->npibonds) {
-		/*  Renumber the pi-atom constructs  */
-		for (i = 0; i < mp->npibonds * 4; i++) {
-			j = mp->pibonds[i];
-			if (j >= 0 && j < mp->natoms)
-				mp->pibonds[i] = old2new[j];
-		}
-	}
-#endif  /*  PIATOM  */
-
 	/*  Renumber the atoms  */
 	for (i = 0; i < mp->natoms; i++)
 		memmove(ATOM_AT_INDEX(mp->atoms, old2new[i]), ATOM_AT_INDEX(saveAtoms, i), gSizeOfAtomRecord);
 	retval = 0;
 	
-#if PIATOM
-	MoleculeInvalidatePiConnectionTable(mp);
-#endif
 	MoleculeIncrementModifyCount(mp);
 	mp->needsMDRebuild = 1;
 
@@ -9832,12 +10082,12 @@ MoleculeCalculateCellFromAxes(XtalCell *cp, int calc_abc)
 		return -1;  /*  Non-regular transform  */
 
 	/*  Calculate the reciprocal cell parameters  */
-	cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[1] * cp->rtr[1] + cp->rtr[2] * cp->rtr[2]);
-	cp->rcell[1] = sqrt(cp->rtr[3] * cp->rtr[3] + cp->rtr[4] * cp->rtr[4] + cp->rtr[5] * cp->rtr[5]);
-	cp->rcell[2] = sqrt(cp->rtr[6] * cp->rtr[6] + cp->rtr[7] * cp->rtr[7] + cp->rtr[8] * cp->rtr[8]);
-	cp->rcell[3] = acos((cp->rtr[3] * cp->rtr[6] + cp->rtr[4] * cp->rtr[7] + cp->rtr[5] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
-	cp->rcell[4] = acos((cp->rtr[6] * cp->rtr[0] + cp->rtr[7] * cp->rtr[1] + cp->rtr[8] * cp->rtr[2]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
-	cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[3] + cp->rtr[1] * cp->rtr[4] + cp->rtr[2] * cp->rtr[5]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
+	cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[3] * cp->rtr[3] + cp->rtr[6] * cp->rtr[6]);
+	cp->rcell[1] = sqrt(cp->rtr[1] * cp->rtr[1] + cp->rtr[4] * cp->rtr[4] + cp->rtr[7] * cp->rtr[7]);
+	cp->rcell[2] = sqrt(cp->rtr[2] * cp->rtr[2] + cp->rtr[5] * cp->rtr[5] + cp->rtr[8] * cp->rtr[8]);
+	cp->rcell[3] = acos((cp->rtr[1] * cp->rtr[2] + cp->rtr[4] * cp->rtr[5] + cp->rtr[7] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
+	cp->rcell[4] = acos((cp->rtr[2] * cp->rtr[0] + cp->rtr[5] * cp->rtr[3] + cp->rtr[8] * cp->rtr[6]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
+	cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[1] + cp->rtr[3] * cp->rtr[4] + cp->rtr[6] * cp->rtr[7]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
 	
 	if (calc_abc) {
 		/*  Calculate a, b, c, alpha, beta, gamma  */
@@ -9848,7 +10098,6 @@ MoleculeCalculateCellFromAxes(XtalCell *cp, int calc_abc)
 		cp->cell[4] = acos((cp->tr[6] * cp->tr[0] + cp->tr[7] * cp->tr[1] + cp->tr[8] * cp->tr[2]) / (cp->cell[2] * cp->cell[0])) * kRad2Deg;
 		cp->cell[5] = acos((cp->tr[0] * cp->tr[3] + cp->tr[1] * cp->tr[4] + cp->tr[2] * cp->tr[5]) / (cp->cell[0] * cp->cell[1])) * kRad2Deg;
 	}
-	
 	return 0;
 }
 
@@ -10036,6 +10285,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;
@@ -10188,16 +10438,6 @@ sMoleculeFragmentSub(Molecule *mp, int idx, IntGroup *result, IntGroup *exatoms)
 			sMoleculeFragmentSub(mp, idx2, result, exatoms);
 		}
 	}
-#if PIATOM
-	if (mp->piconnects != NULL) {
-		for (i = mp->piconnects[idx]; i < mp->piconnects[idx + 1]; i++) {
-			idx2 = mp->piconnects[i];
-			if (IntGroupLookup(result, idx2, NULL))
-				continue;
-			sMoleculeFragmentSub(mp, idx2, result, exatoms);
-		}
-	}
-#endif
 }
 
 /*  The molecular fragment (= interconnected atoms) containing the atom n1 and
@@ -10209,9 +10449,6 @@ MoleculeFragmentExcludingAtomGroup(Molecule *mp, int n1, IntGroup *exatoms)
 	if (mp == NULL || mp->natoms == 0 || n1 < 0 || n1 >= mp->natoms)
 		return NULL;
 	result = IntGroupNew();
-#if PIATOM
-	MoleculeValidatePiConnectionTable(mp);
-#endif
 	sMoleculeFragmentSub(mp, n1, result, exatoms);
 	return result;
 }
@@ -10229,9 +10466,6 @@ MoleculeFragmentExcludingAtoms(Molecule *mp, int n1, int argc, int *argv)
 	for (i = 0; i < argc; i++)
 		IntGroupAdd(exatoms, argv[i], 1);
 	result = IntGroupNew();
-#if PIATOM
-	MoleculeValidatePiConnectionTable(mp);
-#endif
 	sMoleculeFragmentSub(mp, n1, result, exatoms);
 	IntGroupRelease(exatoms);
 	return result;
@@ -10249,9 +10483,6 @@ MoleculeFragmentWithAtomGroups(Molecule *mp, IntGroup *inatoms, IntGroup *exatom
 		return NULL;
 	IntGroupIteratorInit(inatoms, &iter);
 	result = IntGroupNew();
-#if PIATOM
-	MoleculeValidatePiConnectionTable(mp);
-#endif
 	while ((i = IntGroupIteratorNext(&iter)) >= 0) {
 		sMoleculeFragmentSub(mp, i, result, exatoms);
 	}
@@ -10271,9 +10502,6 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
 	Atom *ap;
 	if (mp == NULL || mp->natoms == 0 || group == NULL)
 		return 0;  /*  Invalid arguments  */
-#if PIATOM
-	MoleculeValidatePiConnectionTable(mp);
-#endif
 	bond_count = 0;
 	for (i = 0; (i1 = IntGroupGetStartPoint(group, i)) >= 0; i++) {
 		i2 = IntGroupGetEndPoint(group, i);
@@ -10305,20 +10533,6 @@ MoleculeIsFragmentDetachable(Molecule *mp, IntGroup *group, int *n1, int *n2)
 					}
 				}					
 			}
-#if PIATOM
-			if (mp->piconnects != NULL) {
-				for (k = mp->piconnects[j]; k < mp->piconnects[j + 1]; k++) {
-					k2 = mp->piconnects[k];
-					if (IntGroupLookup(group, k2, NULL) == 0) {
-						bond_count++;
-						nval1 = j;
-						nval2 = k2;
-						if (bond_count > 1)
-							return 0;  /*  Too many bonds  */
-					}
-				}
-			}
-#endif
 		}
 	}
 	if (bond_count == 1) {
@@ -10401,6 +10615,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;
 
@@ -10424,17 +10641,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;
@@ -10448,6 +10662,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,
@@ -10457,23 +10683,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);
@@ -10490,23 +10721,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;
@@ -10522,6 +10767,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)
@@ -10587,38 +10833,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);
@@ -10627,18 +10878,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)
@@ -10650,19 +10903,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;
@@ -10696,7 +10959,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;
@@ -10713,7 +10976,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);
@@ -10738,27 +11001,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 ======
 
-#if !defined(PIATOM)
+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
@@ -10892,134 +11333,27 @@ MoleculeSetPiAnchorList(Molecule *mol, Int idx, Int nentries, Int *entries, Doub
 	return 0;
 }
 
-#endif
-
-#if PIATOM
-int
-MoleculeCalculatePiAtomPosition(Molecule *mol, int idx)
-{
-	PiAtom *pp;
-	Int i, *cp;
-	if (mol == NULL || idx < 0 || idx >= mol->npiatoms)
-		return -1;
-	pp = mol->piatoms + idx;
-	cp = AtomConnectData(&pp->connect);
-	VecZero(pp->r);
-	for (i = pp->connect.count - 1; i >= 0; i--) {
-		Vector rr = ATOM_AT_INDEX(mol->atoms, cp[i])->r;
-		VecScaleInc(pp->r, rr, pp->coeffs[i]);
-	}
-	return idx;
-}
-
-int
-MoleculeValidatePiConnectionTable(Molecule *mol)
-{
-	Int pass, i, j, k, m;
-	
-	if (mol == NULL || mol->pibonds == NULL)
-		return -1; /*  No need to process */
-	if (mol->piconnects != NULL)
-		return 0;  /*  Already valid  */
-
-	/*  Allocate index table  */
-	NewArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), mol->natoms + 1);
-	memset(mol->piconnects, 0, sizeof(Int) * (mol->natoms + 1));
-	
-	/*  Pass 1: count connections for each atom  */
-	/*  Pass 2: store connection info  */
-	for (pass = 0; pass < 2; pass++) {
-		Int n[2], *ip[2], c[2];
-		for (i = 0; i < mol->npibonds; i++) {
-			AtomConnect *ac;
-			if (mol->pibonds[i * 4 + 2] >= 0)
-				continue;  /*  Skip angle or dihedral entries  */
-			for (j = 0; j < 2; j++) {
-				n[j] = mol->pibonds[i * 4 + j];
-				if (n[j] >= ATOMS_MAX_NUMBER) {
-					ac = &(mol->piatoms[n[j] - ATOMS_MAX_NUMBER].connect);
-					ip[j] = AtomConnectData(ac);
-					c[j] = ac->count;
-				} else if (n[j] >= 0 && n[j] < mol->natoms) {
-					ip[j] = &n[j];
-					c[j] = 1;
-				} else break;
-			}
-			if (j < 2)
-				continue;  /*  Ignore the invalid entry  */
-			for (j = 0; j < c[0]; j++) {
-				for (k = 0; k < c[1]; k++) {
-					Int a1 = ip[0][j];
-					Int a2 = ip[1][k];
-					if (pass == 0) {
-						/*  Count  */
-						mol->piconnects[a1]++;
-						mol->piconnects[a2]++;
-					} else {
-						/*  Store the entry (at the first empty slot)  */
-						for (m = mol->piconnects[a1]; m < mol->piconnects[a1 + 1]; m++) {
-							if (mol->piconnects[m] == -1) {
-								mol->piconnects[m] = a2;
-								break;
-							}
-						}
-						for (m = mol->piconnects[a2]; m < mol->piconnects[a2 + 1]; m++) {
-							if (mol->piconnects[m] == -1) {
-								mol->piconnects[m] = a1;
-								break;
-							}
-						}
-					}
-				}
-			}
-		}
-		if (pass == 0) {
-			/*  Expand the table, and store the position numbers  */
-			m = mol->natoms + 1;
-			for (i = 0; i <= mol->natoms; i++) {
-				j = mol->piconnects[i];
-				mol->piconnects[i] = m;
-				m += j;
-			}
-			AssignArray(&mol->piconnects, &mol->npiconnects, sizeof(Int), m - 1, NULL);
-			for (j = mol->natoms + 1; j < m; j++)
-				mol->piconnects[j] = -1;
-		}
-	}
-	return (mol->npiconnects - mol->natoms - 1);  /*  Returns the number of entries  */
-}
-
-void
-MoleculeInvalidatePiConnectionTable(Molecule *mol)
-{
-	if (mol == NULL)
-		return;
-	if (mol->piconnects != NULL) {
-		free(mol->piconnects);
-		mol->piconnects = NULL;
-	}
-	mol->npiconnects = 0;
-}
-#endif  /*  PIATOM  */
-
 #pragma mark ====== MO calculation ======
 
 /*  Calculate an MO value for a single point.  */
-/*  Index is the MO number (1-based)  */
+/*  Index is the MO number (1-based); 0 denotes "arbitrary vector"  */
 /*  tmp is an array of (natoms * 4) atoms, and used to store dr and |dr|^2 for each atom.  */
 static Double
-sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
+sCalcMOPoint(Molecule *mp, const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
 {
 	ShellInfo *sp;
 	PrimInfo *pp;
 	Double val, tval, *cnp, *tmpp, *mobasep, *mop;
 	Int i, j;
 	/*  Cache dr and |dr|^2  */
-	for (i = 0; i < bset->natoms; i++) {
-		Vector r = bset->pos[i];
-		tmp[i * 4] = r.x = vp->x - r.x;
-		tmp[i * 4 + 1] = r.y = vp->y - r.y;
-		tmp[i * 4 + 2] = r.z = vp->z - r.z;
+	if (index == 0)
+		index = bset->nmos + 1;
+	for (i = 0; i < mp->natoms; i++) {
+		Vector r;
+		r = ATOM_AT_INDEX(mp->atoms, i)->r;
+		tmp[i * 4] = r.x = (vp->x - r.x) * kAngstrom2Bohr;
+		tmp[i * 4 + 1] = r.y = (vp->y - r.y) * kAngstrom2Bohr;
+		tmp[i * 4 + 2] = r.z = (vp->z - r.z) * kAngstrom2Bohr;
 		tmp[i * 4 + 3] = r.x * r.x + r.y * r.y + r.z * r.z;
 	}
 	/*  Iterate over all shells  */
@@ -11028,6 +11362,8 @@ sCalcMOPoint(const BasisSet *bset, Int index, const Vector *vp, Double *tmp)
 	for (i = 0, sp = bset->shells; i < bset->nshells; i++, sp++) {
 		pp = bset->priminfos + sp->p_idx;
 		cnp = bset->cns + sp->cn_idx;
+		if (sp->a_idx >= mp->natoms)
+			return 0.0; /*  This may happen when molecule is edited after setting up MO info  */
 		tmpp = tmp + sp->a_idx * 4;
 		mop = mobasep + sp->m_idx;
 		switch (sp->sym) {
@@ -11116,13 +11452,14 @@ 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  */
 		}
 	}
 	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)
 {
@@ -11135,6 +11472,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;
@@ -11154,7 +11494,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);
@@ -11168,7 +11508,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;
@@ -11187,35 +11527,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;
@@ -11292,17 +11635,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);
 	
@@ -11310,7 +11657,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");
 			}
@@ -11319,3 +11679,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);
+}