Rebase/softcore v2 onto main#110
Open
David-Araripe wants to merge 124 commits into
Open
Conversation
…optimal bond geometries from auto-added hydrogens
…ns, solute. Now, setting the FEP run to 2fs will apply shake to all hydrogens in the simulation. Solvent SHAKE is set to off in case solvent contains other solvents like POPC. shake_hydrogens already handles shake settings for HOH molecules.
- add generate_golden_restraints.py to generate restraints compatible with QligFEP v2.1.0
…ame to avoid conflict with python package name (avoid conflicts when importing from pytest directory)
… setter is doing its job correctly
…oduction runs. This commit includes the additon of a "templates" subpackage that will handle the input files instead of the old INPUTS/. An example on how to display the rendered MD files was added to templates.__init__
…ized atom type, parse the whole protein file and tell all the problematic atom/residue names
…igfep. Qprep input sare now passed as positional arguments allowing us to use subprocess for command execution. qprep positional arguments are supported commit 8d67120
1. use the qprep & qfep templates from the templates/ directory
2. dynamicall remove out-of-sphere molecules that aren't part of the system (avoid unecessary calculation of bonded energies with the MD engine)
Also added tests for this fragment filtering step along with two protein structures that serve as test cases
…n pair list construction Build compressed sparse row (CSR) lookup tables from listexlong and list14long at startup. Each atom gets a sorted list of its long-range exclusion/1-4 partners, enabling binary search instead of scanning the full list for every atom pair. Replaces linear scans in all 8 nbpplis2*/nbpplist* routines plus nbpp_count and the monitor group 1-4 check. Benchmarks (5000 steps, eq1): - CDK2 protein: 109s -> 96s (1.13x speedup) - P2Y1 GPCR membrane (LRF): 488s -> 179s (2.73x speedup) - Tyk2 water: 70s -> 70s (no regression) Bit-for-bit identical energies on all three systems.
…ing a single approach, the method now attempts to generate the conformer in a tirered approach: 1. ConstrainedEmbed with element-strict MCS 2. ConstrainedEmbed with topology-only MCS (only if it differs from #1) 3. coordMap embed with topology-only MCS positions 4. Unconstrained embed + centroid translation
…x alingment issues and repair the network and edge topology
…routine (out-of-sphere charges shouldn't be present)
…ffset. It applies Beutler 1994 Eq. 7's (1-lambda_s)^2 modulation to the offset itself per state, so that: lambda_s = 1: offset = 0 -> standard LJ + Coulomb lambda_s = 0: offset = full -> bounded V at r=0 for cross-state evaluation intermediate: offset = full * (1 - lambda_s)^2 Scaling the offset (post-quadratic), not the user's max-V target, preserves endpoint behavior - target-scaling would propagate through the quadratic in populate_sc_lookup and leave a non-vanishing offset at lambda=1. The Coulomb term uses the s=p=6 case of Beutler's generalized Eq. 9 so LJ and Coulomb can share a single offset table; this is documented in the softcore_scale_beutler docstring. QligFEP now writes symmetric (20, 20) alphas when softcore_method=beutler_coul so the modulation can protect both endpoints; the legacy asymmetric (0, 20) pattern is preserved for back compatibility. Refactor: the sc_lookup fill loop is extracted from qatom_load_fep into populate_sc_lookup so softcore_scale_beutler can reuse it.
…1.1-1.4 The previous SC_GAPSYS implementation had three correctness bugs against the source manuscript (Gapsys, Seeliger, de Groot, JCTC 2012, 8, 2373) and its Supporting Information §S1: 1. Lambda inversion in softcore_init_gapsys. The paper formulas assume lambda_state increases as the state disappears. Q's convention is the opposite (EQ(istate)%lambda = 1 means coupled/real, = 0 means decoupled/ghost), so r_sc must scale with (1 - lambda_state)^(1/6), not lambda_state^(1/6). The previous code applied no softcore at the ghost endpoint exactly when it was needed -- the root cause of intra-ligand Q-Q endpoint overflow under Gapsys. 2. Endpoint-skip cycle (if lambda_st < 1e-10 cycle) compounded #1 by leaving the lookup arrays at zero for the ghost state. Removed; the corrected formula handles the endpoint naturally. 3. Runtime V(r), F(r) inside r < r_sc did not match SI Eqs. 1.1-1.4. The previous Q form was continuous but qualitatively wrong (F(0) = 0, force increasing with r). Paper requires F(0) = 168 C12/r_sc^13 - 48 C6/r_sc^7 and force monotonically decreasing with r (paper §2.3). Replaced with the SI's exact a*r^2 - b*r + c form for LJ and Coulomb via two new pure helpers gapsys_eval_lj and gapsys_eval_q in qatom.f90, called from each of the ~15 Gapsys runtime sites in md.f90 and md_test.f90. Also addressed: - Exposed sigma_Q as a tunable parameter (FEP keyword gapsys_sigma_q, default 1.0 per paper). The previously-existing gapsys_sigma_lj is unrelated to the paper's sigma_Q -- it is the LJ sigma fallback for dummy atoms with C6 = C12 = 0. - Q-protein and Q-water Coulomb radius now uses the full paper Eq. 5 factor (1 + sigma_Q * |q_i*q_j|) at runtime where q_j is in scope. Previously precomputation used |q_i| only because partner charges are per-atom, not per-type. Restructured: precompute alpha_Q*(1-lambda)^(1/6) per (q-atom, state) as gapsys_q_lambda_factor (2D), multiply by the partner-charge factor at each nonbonded runtime site. Old gapsys_q_rsc (3D) removed. - R*-to-sigma conversion under arithmetic mixing. AMBER14sb.prm and CHARMM36.prm store R* (R_min/2), not sigma/2. Q's standard LJ formula absorbs the R_min-vs-sigma factor cleanly because it consumes R_min^12 and 2*R_min^6 directly, but Gapsys's r_sc formula is in terms of sigma (zero-crossing radius). The arithmetic branch of softcore_init_gapsys now divides the summed R* values by 2^(1/6) to obtain sigma_pair. Without this, r_sc was ~12% too large under AMBER and CHARMM. The geometric branch is unaffected -- (C12/C6)^(1/6) extracts sigma_pair directly from physical C12/C6. - MPI broadcast of all four Gapsys scalars (alpha_lj, alpha_q, sigma_lj, sigma_q) plus the new gapsys_q_lambda_factor 2D array. The previous block broadcast only the 3D lookup tables, not the scalars, so non-master ranks held stale defaults. - qligfep.py comment updated: the asymmetric (0,20)/(20,0) [softcore] alpha pattern for Gapsys is now redundant rather than load-bearing. The (1-lambda)^(1/6) formula already produces r_sc = 0 at the coupled endpoint, so the sc_lookup gate is a no-op there. Pattern kept for backward compatibility. Out of scope (deferred to follow-up): - qfep TI estimator. SI Eqs. 1.5-1.6 give the full dH/dlambda including the chain-rule term from the lambda-dependence of r_sc(lambda). qfep currently approximates dH/dlambda from per-state V_1, V_2 only, dropping that contribution. BAR and OS use per-state V directly and are now fully paper-faithful; TI awaits a follow-up. Files changed: src/q6/qatom.f90 new gapsys_sigma_q param + parsing; rewritten softcore_init_gapsys with (1-lambda)^(1/6) convention, R*-to-sigma conversion, and gapsys_q_lambda_factor replacing gapsys_q_rsc; new pure helpers gapsys_eval_lj and gapsys_eval_q implementing SI Eqs. 1.1-1.4. src/q6/md.f90 all Gapsys runtime sites call the new helpers and apply the partner-charge factor at runtime; MPI block updated. src/q6/md_test.f90 mirrored md.f90 changes. src/QligFEP/qligfep.py comment-only update on the [softcore] writer.
…al water and ions on the water leg to match the protein leg total charge; still to be tested)
…lambda 0.988 and 0.013 to be written, for example
…Changes include:
pdb2amber / atom renaming:
- drop the unconditional H1->H / H2->H mappings that corrupted real NH3+ N-termini
(ARG/GLU/GLN/THR/LEU/GLY); add context-aware lone-backbone-H handling instead
- NME cap methyl: map HH31/HH32/HH33 -> H1/H2/H3
- ASN sidechain HD1/HD2 -> HD21/HD22
- generalize ion renaming and make the lookup case-insensitive; add Ca -> CAL
(factor_xa ships calcium as residue "Ca" / atom "CA1")
Force field (AMBER14sb):
- add {SEP} (phosphoserine) library entry, generated with qtools/QLib from
AmberTools phosaa14SB.lib (same source as TPO); net charge -2.00. All atom types
and bonded params were already present from the TPO import.
- add {CAL} ion + CAL atom type, Li & Merz 2013 CM set for TIP3P
(frcmod.ions234lm_126_tip3p, R*=1.649, eps=0.10592870) -- the set MAG/ZIN use.
FF/README.md:
- fix the atom_types column-5 rule: it is the 1-4 epsilon (eps/2), not sigma/2.
The [atom_types] column 5 is the 1-4 Lennard-Jones well depth (epsilon2&3): qprep parses it into bvdw(3) (prep.f90), md sqrt's it (md.f90), and every 1-4 pair uses it via LJcod=3 (calc_nb.f90, md.f90), giving eps_pair_14 = sqrt(col5_i*col5_j) in the arithmetic LJ formula. AMBER scales 1-4 vdW by 1/SCNB = 1/2 at the same R*, so column 5 must be eps/2. 31 hand-added types stored R*/2 there instead (inheriting the wrong rule the FF README used to state): 8 phospho-AA (OP OQ OR OT OV OX OZ CG), 20 POPC Lipid21 (lowercase), 3 ions (CHL MAG ZIN). Their normal LJ (col1/col3) matches the AMBER source exactly (frcmod.phosaa14SB, parm10.dat) -- only col5 was mistyped. This made phospho-AA and POPC intra-molecular 1-4 vdW several-fold too strong; the ions were unaffected (monatomic, no 1-4 pairs).
The Neutralizer identified residues to rewrite by (chain_id, residue_seq_number) only. Two residues can share that pair via an insertion code -- classically a C-terminal NME cap numbered 167A on top of residue 167. Neutralizing residue 167 then also rewrote the cap (renamed it / stripped atoms), so qprep could not build it and excluded the whole protein (complexnotexcluded.pdb came out water-only). Add insertion_code to the match mask in all three rewrite paths (general neutralization, N-terminal, C-terminal) and carry it into the residue-info dict. Restores preparation for the 14 IndustryBenchmarks targets with capped residues sharing chain+seq via insertion codes (back to 60/60). Test: TestNeutralizerInsertionCode in test_qprep_cli.py.
…rameters. Golden files are changed to state "polarized", the actual spelling used by src/q6/md.f90. The word "polarised" was silently ignored, setting the parameters to default regardless of input value
…ng in the same exact coordinate is now checked in a while loop instead of a single increment of 0.001 to the coordinate
…f the legacy `replacements` dictionary-based approach
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Rebase softcore-methods-v2 onto main
Brings the
feature/softcore-methods-v2line of work up to date on top ofmain(124 commits). This is a broad update spanning the MD engine, the QligFEP setup layer, perturbation-network generation, system preparation, and the test suite. Highlights below, grouped by theme.Soft-core methods
Adds a selectable soft-core scheme via
--softcore-method {standard, beutler_coul, gapsys}(engine constantsSC_STANDARD,SC_BEUTLER_COUL,SC_GAPSYS), threaded from the CLI through the FEP-file writer into all 15 nonbonded routines of the engine:standard— soft-core on LJ only (default; prior behaviour).beutler_coul— extends soft-core to Coulomb via a modified effective distance, with the Beutler-1994(1-λ)²offset scaling.gapsys— force-based linearization below a critical radius for both LJ and Coulomb, matched to Gapsys–Seeliger–de Groot (2012) Eqs. 1.1–1.4.Mirrored into
md_test.f90, initialized at FEP setup (softcore_init_gapsys), and covered by new tests.Important
All three methods are in place and selectable, but the way Q builds dual-topology perturbations means none of them are actually triggered at runtime. The code paths exist, are threaded through the engine, and are exercised by tests — but Q's per-state Hamiltonian mixing with the dual-topology DUM convention leaves the soft-core offset structurally dormant: it evaluates to zero on the DUM (vanishing) side, and QligFEP's asymmetric
[softcore]block disables it on the real side, so the offset never modifies the forces seen by the integrator. The three methods therefore produce statistically equivalent trajectories for routine congeneric FEP.Perturbation-network generation
qkonnektorCLI (konnektor-based networks) alongside the existingqlomap, with aRestraintScorerthat combines LOMAP + RestraintSetter for better edge quality.--central_ligand(star networks), redundant-MST and n-edges network types,--self-solve(repair alignment/topology),--exp-key(standardized experimental-value naming),--separate_charges.GlobalLigandAlignerrenamed tofkcombuLigandAligner(clearer), withscaffold_lock, atom/bond-type options, and anoptions()inspector. TieredBiasedConformerGeneratorwith topology-only MCS default.qlomap/qkonnektorto match Q's thermodynamic cycle.Template system refactor (partial migration)
Moves input-file generation off the legacy
replacements-dictionary approach toward render-based templates intemplates/:qprep, andqfepinputs now render fromtemplates/with structured string formatting.replacementspath still exists; this PR moves the codebase toward the template approach rather than removing the old one. Golden-file "identical generation" tests lock the rendered output.Charge-changing perturbations
System preparation robustness (qprep / pdb2amber)
pdb2amber+ AMBER14sb hardened across many targets: proline handling,ATOM/HETATM-only parsing,TER-sensitivity.qprep_cliand QligFEP; hydrogen-placement optimization; POPC lipid charge-groups for sphere exclusion.MD engine — minimization, performance, MPI
[MD]parameters;eq1default updated to the new protocol.Notable bugfixes
0.988/0.013instead of the intended values.polarizedvspolarised: golden files corrected to the spelling the engine actually reads (polarisedwas silently ignored, defaulting the parameters).qprep.f90positional-argument bug; variousUnboundLocalError/NameError/duplicate-logging fixes.Charges, dependencies, tests
--nagl-modelpassthrough.Notes for reviewers
qdyn_neqengine intoqdyn/qdynpis a separate branch (refactor/consolidate-neq-engine) that will be rebased onto main and opened as its own PR after this one lands.