add tether force, group_constraint_by_position

functionTerms/constraints.py

@@ -17,3 +17,46 @@ def group_constraint_by_distance(oa, d0=0*angstrom, group1=[oa.ca[0], oa.ca[1]],
     constraint.addBond([0, 1])
     constraint.setForceGroup(forceGroup)
     return constraint
+
+from simtk.unit import dalton
+def group_constraint_by_position(oa, k=1*kilocalorie_per_mole, x0=10, y0=10, z0=10, appliedToResidues=-1, forceGroup=3):
+    # x0, y0, z0 is in unit of nm.
+    # appliedToResidues can be a list of residue index. for example appliedToResidues=[0, 1], to tether the first two residues.
+    # 1 Kcal = 4.184 kJ strength by overall scaling
+    k = k.value_in_unit(kilojoule_per_mole)   # convert to kilojoule_per_mole, openMM default uses kilojoule_per_mole as energy.
+    k_constraint = k * oa.k_awsem
+    sum_of_x_coord = CustomExternalForce(f"x*mass")
+    sum_of_y_coord = CustomExternalForce(f"y*mass")
+    sum_of_z_coord = CustomExternalForce(f"z*mass")
+
+    sum_of_x_coord.addPerParticleParameter("mass")
+    sum_of_y_coord.addPerParticleParameter("mass")
+    sum_of_z_coord.addPerParticleParameter("mass")
+
+    # print("index for CAs", oa.ca)
+    print(f"mass can be retrieved as ", oa.system.getParticleMass(oa.ca[0]))
+    total_mass = 0.0
+    for i in range(oa.natoms):
+        if appliedToResidues == -1:
+            mass = oa.system.getParticleMass(i).value_in_unit(dalton)
+            sum_of_x_coord.addParticle(i, [mass])
+            sum_of_y_coord.addParticle(i, [mass])
+            sum_of_z_coord.addParticle(i, [mass])
+            total_mass += mass
+        elif oa.resi[i] in appliedToResidues:
+            mass = oa.system.getParticleMass(i).value_in_unit(dalton)
+            sum_of_x_coord.addParticle(i, [mass])
+            sum_of_y_coord.addParticle(i, [mass])
+            sum_of_z_coord.addParticle(i, [mass])
+            total_mass += mass
+        # if oa.resi[i] == appliedToResidue:
+        #     pulling.addParticle(i)
+        # print(oa.resi[i] , oa.seq[oa.resi[i]])
+    print(f"total_mass = {total_mass}")
+    harmonic = CustomCVForce(f"{k_constraint}*((sum_x/{total_mass}-{x0})^2+(sum_y/{total_mass}-{y0})^2+(sum_z/{total_mass}-{z0})^2)")
+    harmonic.addCollectiveVariable("sum_x", sum_of_x_coord)
+    harmonic.addCollectiveVariable("sum_y", sum_of_y_coord)
+    harmonic.addCollectiveVariable("sum_z", sum_of_z_coord)
+    harmonic.setForceGroup(forceGroup)
+    return harmonic
+

functionTerms/qBiasTerms.py

@@ -149,11 +149,12 @@ def partial_q_value(oa, reference_pdb_file, min_seq_sep=3, a=0.1, startResidueIn
     qvalue.setForceGroup(forceGroup)
     return qvalue
 
-def qbias_term(oa, q0, reference_pdb_file, reference_chain_name, k_qbias=100*kilocalorie_per_mole, qbias_min_seq_sep=3, qbias_max_seq_sep=np.inf, qbias_contact_threshold=0.8*nanometers):
+def qbias_term(oa, q0, reference_pdb_file, reference_chain_name, k_qbias=100*kilocalorie_per_mole, qbias_min_seq_sep=3, qbias_max_seq_sep=np.inf, qbias_contact_threshold=0.8*nanometers, forceGroup=4):
     k_qbias = k_qbias.value_in_unit(kilojoule_per_mole)   # convert to kilojoule_per_mole, openMM default uses kilojoule_per_mole as energy.
-    qbias = CustomCVForce("0.5*k_qbias*(q-q0)^2")
+    qbias = CustomCVForce(f"0.5*{k_qbias}*(q-{q0})^2")
     q = q_value(oa, reference_pdb_file, reference_chain_name, min_seq_sep=qbias_min_seq_sep, max_seq_sep=qbias_max_seq_sep, contact_threshold=qbias_contact_threshold)
     qbias.addCollectiveVariable("q", q)
-    qbias.addGlobalParameter("k_qbias", k_qbias)
-    qbias.addGlobalParameter("q0", q0)
+    # qbias.addGlobalParameter("k_qbias", k_qbias)
+    # qbias.addGlobalParameter("q0", q0)
+    qvalue.setForceGroup(forceGroup)
     return qbias

mm_analysis.py

@@ -140,6 +140,7 @@
                     # , "Q_wat":4, "Q_mem":5, "Debye_huckel":30
                    }
 print("Please ensure the forceGroupTable in mm_analysis is set up correctly if you are adding new energy terms.")
+print("Also, please notice that the total energy include all the terms with group index range from 11 to 32.")
 # forceGroupTable = {"Con":11, "Chain":12, "Chi":13, "Excluded":14, "Rama":15, "Direct":16,
 #                     "Burial":17, "Mediated":18, "Contact":18, "Fragment":19, "Membrane":20, "ER":21,"TBM_Q":22, "beta_1":23, "beta_2":24,"beta_3":25,"pap":26, "Total":list(range(11, 32)),
 #                     "Water":[16, 18], "Beta":[23, 24, 25], "Pap":26, "Rg_Bias":27, "Helical":28, "Pulling":29, "Q":1, "Rg":2, "Qc":3, "Q_wat":4, "Q_mem":5}