Merge pull request #1454 from virtualcell/dan-ss-rates-sanity

Langevin input file gets wrong particle location name

Author: danv61

vcell-client/src/main/java/cbit/vcell/client/desktop/biomodel/ReactionRulePropertiesTableModel.java

@@ -241,6 +241,9 @@ public boolean isCellEditable(int rowIndex, int columnIndex) {
 	 *   and the property that has changed.
 	 */
 public void propertyChange(java.beans.PropertyChangeEvent evt) {
+	// there is also an "entityChange" propertyChange event generated from ReactionRulePropertiesTreeModel (which is deprecated)
+	// when a rule is selected in the upper panel (BioModelEditorReactionTableModel)
+	// TODO: make sure we ignore it everywhere, should be phased out
 	if (evt.getSource() == this && evt.getPropertyName().equals("reactionRule")) {
 		ReactionRule oldValue = (ReactionRule)evt.getOldValue();
 		if (oldValue != null){

vcell-client/src/main/java/cbit/vcell/mapping/gui/MolecularStructuresPropertiesPanel.java

@@ -232,9 +232,23 @@ public void mouseClicked(MouseEvent e) {
             gbc.anchor = GridBagConstraints.WEST;
             optionsPanel.add(getZoomSmallerButton(), gbc);
 
+//            gbc = new GridBagConstraints();
+//            gbc.gridx = 0;
+//            gbc.gridy = 2;
+//            gbc.anchor = GridBagConstraints.WEST;
+//            gbc.insets = new Insets(4, 4, 4, 10);
+//            optionsPanel.add(new JLabel("Reaction Radius"), gbc);
+//
+//            gbc = new GridBagConstraints();
+//            gbc.gridx = 1;
+//            gbc.gridy = 2;
+//            gbc.anchor = GridBagConstraints.WEST;
+//            gbc.insets = new Insets(4, 4, 4, 10);
+//            optionsPanel.add(new JLabel("2 nm"), gbc);
+
             gbc = new GridBagConstraints();
             gbc.gridx = 0;
-            gbc.gridy = 2;
+            gbc.gridy = 2;          // make this 3 if we show reaction radius here
             gbc.weightx = 1;
             gbc.weighty = 1;		// fake cell used for filling all the vertical empty space
             gbc.anchor = GridBagConstraints.WEST;

vcell-client/src/main/java/org/vcell/model/springsalad/gui/MolecularStructuresPanel.java

@@ -49,6 +49,7 @@
 import java.awt.event.FocusEvent;
 import java.awt.event.FocusListener;
 import java.beans.PropertyChangeListener;
+import java.text.DecimalFormat;
 import java.util.ArrayList;
 import java.util.Map;
 
@@ -543,7 +544,12 @@ public Component getTableCellRendererComponent(JTable table, Object value, boole
 								} else {
 									setText(value + "");		// if it's too busy, just show the numbers
 								}
-								setToolTipText(value + " [nm]");	// we always show the units in the tooltip
+								double minRadius = Math.max(0.5+dvalue,  1.5*dvalue);
+								Double R = Math.min(minRadius, dvalue+2);
+								double truncatedR = Math.floor(R * 100) / 100.0;
+								DecimalFormat decimalFormat = new DecimalFormat("#.##");
+								String formattedR = decimalFormat.format(truncatedR);
+								setToolTipText(value + " [nm], (Reaction Radius: " + formattedR + " [nm])");	// we always show the units in the tooltip
 								break;
 							case COLUMN_DIFFUSION:
 								if(cellWidth > 70) {

vcell-core/src/main/java/cbit/vcell/mapping/ReactionRuleSpec.java

@@ -1201,12 +1201,14 @@ public void gatherIssues(IssueContext issueContext, List<Issue> issueList, React
 				if(scCandidate == null || scCandidate.getSpeciesPattern() == null) {
 					continue;
 				}
+				// we may have just one spCandidate wif the reaction is A + A -> A.A
 				SpeciesPattern spCandidate = scCandidate.getSpeciesPattern();
 				MolecularTypePattern mtpCandidate = spCandidate.getMolecularTypePatterns().get(0);
 				MolecularType mtCandidate = mtpCandidate.getMolecularType();
 				if(mtOursOne == mtCandidate) {
 					siteAttributesMapOne = scs.getSiteAttributesMap();
-				} else if(mtOursTwo == mtCandidate) {
+				}
+				if(mtOursTwo == mtCandidate) {
 					siteAttributesMapTwo = scs.getSiteAttributesMap();
 				}
 			}
@@ -1225,13 +1227,14 @@ public void gatherIssues(IssueContext issueContext, List<Issue> issueList, React
 				}
 			}
 			//
-			// TODO: check logic here, sometimes siteAttributesMapTwo is null (issue #977)
-			// temp fix, we check for non-null siteAttributesMapTwo
-			// hard to catch, it may be that the reactant or the product end up in an inconsistent state
-			// after using the reaction visual editor in the physiology (the state of a bonding site is possibly inconsistent)
-			//
-			if(sasOne != null && sasTwo != null && siteAttributesMapTwo != null) {
-				
+			// with previous logic, siteAttributesMapTwo could have been null (issue #977)
+			// happened when the reaction was between 2 molecules of the same type  A + A -> A.A
+			// the temp fixwas:  we check for non-null siteAttributesMapTwo
+			// as from feb 24, 2025 this was fixed and siteAttributesMapTwo should never be null
+			if(siteAttributesMapTwo == null) {
+				throw new RuntimeException("Unexpected null value for siteAttributesMapTwo");
+			}
+			if(sasOne != null && sasTwo != null) {
 				for(Map.Entry<MolecularComponentPattern, SiteAttributesSpec> entry : siteAttributesMapTwo.entrySet()) {
 					MolecularComponentPattern mcpCandidate = entry.getKey();
 					if(MolecularComponentPattern.BondType.None != mcpCandidate.getBondType()) {
@@ -1250,12 +1253,12 @@ public void gatherIssues(IssueContext issueContext, List<Issue> issueList, React
 					issueList.add(new Issue(r, issueContext, IssueCategory.Identifiers, msg, tip, Issue.Severity.ERROR));
 					return;
 				}
-//				if(checkOnRate(sasOne, sasTwo) == false) {	// rate doesn't check as acceptable
-//					String msg = "The simulation Kon is too large (I.e. exceeds the diffusion limited rate) for this reaction.";
-//					String tip = "Please consider reducing Kon or increasing the Radius or D of the participating Site Types.";
-//					issueList.add(new Issue(r, issueContext, IssueCategory.Identifiers, msg, tip, Issue.Severity.ERROR));
-//					return;
-//				}
+				if(checkOnRate(sasOne, sasTwo) == false) {	// rate doesn't check as acceptable
+					String msg = "The forward rate Kf is too large (i.e. exceeds the diffusion limited rate) for this reaction rule.";
+					String tip = "Please consider reducing Kon or increasing the Radius or D of the participating Site Types.";
+					issueList.add(new Issue(r, issueContext, IssueCategory.Identifiers, msg, tip, Issue.Severity.WARNING));
+					return;
+				}
 			}
 			// binding reactions must be reversible
 			if(!reactionRule.isReversible()) {
@@ -1321,6 +1324,11 @@ public void gatherIssues(IssueContext issueContext, List<Issue> issueList, React
 	 */
 	public boolean checkOnRate(SiteAttributesSpec sasOne, SiteAttributesSpec sasTwo) {
 
+		// set of acceptable numbers (marginally) for A + A -> A.A are:
+		// Kon = 40 s-1uM-1
+		// site radius 1nm
+		// site diffusion rate 1 um2s-1
+
 		double R = sasOne.computeReactionRadius() + sasTwo.computeReactionRadius();	// nm
 		double D = sasOne.getDiffusionRate() + sasTwo.getDiffusionRate();
 

vcell-core/src/main/java/org/vcell/solver/langevin/LangevinLngvWriter.java

@@ -5,6 +5,7 @@
 
 import cbit.vcell.geometry.AnalyticSubVolume;
 import cbit.vcell.geometry.SubVolume;
+import cbit.vcell.math.*;
 import cbit.vcell.model.Structure;
 import cbit.vcell.solver.*;
 import org.vcell.util.Pair;
@@ -18,28 +19,8 @@
 import cbit.vcell.mapping.SimulationContext;
 import cbit.vcell.mapping.SpeciesContextSpec;
 import cbit.vcell.mapping.SimulationContext.Application;
-import cbit.vcell.math.Action;
-import cbit.vcell.math.JumpProcessRateDefinition;
-import cbit.vcell.math.LangevinParticleJumpProcess;
-import cbit.vcell.math.LangevinParticleMolecularComponent;
-import cbit.vcell.math.LangevinParticleMolecularType;
-import cbit.vcell.math.MacroscopicRateConstant;
-import cbit.vcell.math.MathDescription;
-import cbit.vcell.math.MathException;
-import cbit.vcell.math.MathUtilities;
-import cbit.vcell.math.ParticleComponentStateDefinition;
-import cbit.vcell.math.ParticleComponentStatePattern;
-import cbit.vcell.math.ParticleJumpProcess;
-import cbit.vcell.math.ParticleMolecularComponent;
-import cbit.vcell.math.ParticleMolecularComponentPattern;
-import cbit.vcell.math.ParticleMolecularType;
-import cbit.vcell.math.ParticleMolecularTypePattern;
-import cbit.vcell.math.ParticleProperties;
 import cbit.vcell.math.ParticleProperties.ParticleInitialCondition;
 import cbit.vcell.math.ParticleProperties.ParticleInitialConditionCount;
-import cbit.vcell.math.ParticleSpeciesPattern;
-import cbit.vcell.math.SubDomain;
-import cbit.vcell.math.Variable;
 import cbit.vcell.mathmodel.MathModel;
 import cbit.vcell.messaging.server.SimulationTask;
 import cbit.vcell.parser.DivideByZeroException;
@@ -910,7 +891,21 @@ private static void writeSpeciesInfo(StringBuilder sb) {
 				throw new RuntimeException("Initial concentration must be a number");
 			}
 
-			sb.append("MOLECULE: \"" + lpmt.getName() + "\" " + subDomain.getName() +
+			String locationName;
+			if(subDomain instanceof MembraneSubDomain) {
+				// we can't use the "Extracellular_Intracellular_membrane" math name of the subdomain,
+				// the Langevin solver has no clue about what that might be
+				// luckily we enforce that a membrane molecule must be anchored to exactly one membrane (named "Membrane")
+				// alternatively we can check the Anchor site and get the membrane name "Membrane" from there
+				if(lpmt.getAnchorList().size() != 1) {
+					throw new RuntimeException("MembraneSubDomain particle must be anchored to one and only one Membrane");
+				}
+				// we make sure elsewhere that it's a membrane
+				locationName = lpmt.getAnchorList().get(0);
+			} else {
+				locationName = subDomain.getName();
+			}
+			sb.append("MOLECULE: \"" + lpmt.getName() + "\" " + locationName +
 					" Number " + scount + 
 					// number of site types and number of sites is the same for the vcell implementation of springsalad
 					" Site_Types " + lpmt.getComponentList().size() + " Total"  + "_Sites " + lpmt.getComponentList().size() +