<commit_before>#!/usr/local/bin/env python

"""
Test YANK using simple models.

DESCRIPTION

This test suite generates a number of simple models to test the 'Yank' facility.

COPYRIGHT

Written by <NAME> <<EMAIL>> while at the University of California Berkeley.

LICENSE

This code is licensed under the latest available version of the GNU General Public License.

"""

#=============================================================================================
# GLOBAL IMPORTS
#=============================================================================================

import os
import sys
import math
import copy
import time
import datetime
from functools import partial

import numpy

from simtk import openmm
from openmmtools import testsystems

from yank import Yank
from yank.repex import ThermodynamicState

import logging
logger = logging.getLogger(__name__)

#=============================================================================================
# MODULE CONSTANTS
#=============================================================================================

from simtk import unit
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA # Boltzmann constant

#=============================================================================================
# MAIN AND TESTS
#=============================================================================================

def test_LennardJonesPair():
    """
    Compute binding free energy of two Lennard-Jones particles and compare to numerical result.

    """

    NSIGMA_MAX = 6.0 # number of standard errors tolerated for success

    # Create Lennard-Jones pair.
    test = testsystems.LennardJonesPair()
    system, positions = test.system, test.positions
    thermodynamic_state = ThermodynamicState(temperature=300.0*unit.kelvin)
    binding_free_energy = test.get_binding_free_energy(thermodynamic_state)
    kT = kB * thermodynamic_state.temperature

    # Create temporary directory for testing.
    import tempfile
    store_dir = tempfile.mkdtemp()
    
    # Initialize YANK object.
    options = dict()
    options['restraint_type'] = 'flat-bottom'
    options['number_of_iterations'] = 50
    options['platform'] = openmm.Platform.getPlatformByName("Reference") # use Reference platform for speed
    options['mc_rotation'] = False
    options['mc_displacement'] = True
    options['mc_displacement_sigma'] = 1.0 * unit.nanometer
    options['verbose'] = False
    options['timestep'] = 5 * unit.femtoseconds
    options['nsteps_per_iteration'] = 20

    # Override box vectors.
    box_edge = 30.0 * unit.angstrom
    a = unit.Quantity((box_edge,        0 * unit.angstrom, 0 * unit.angstrom))
    b = unit.Quantity((0 * unit.angstrom, box_edge,        0 * unit.angstrom))
    c = unit.Quantity((0 * unit.angstrom, 0 * unit.angstrom, box_edge))
    system.setDefaultPeriodicBoxVectors(a, b, c)

    # Alchemical protocol.
    from yank.alchemy import AlchemicalState
    alchemical_states = list()
    lambda_values = [0.0, 0.25, 0.50, 0.75, 1.0]
    for lambda_value in lambda_values:
        alchemical_state = AlchemicalState()
        alchemical_state['lambda_electrostatics'] = lambda_value
        alchemical_state['lambda_sterics'] = lambda_value
        alchemical_states.append(alchemical_state)
    protocols = dict()
    protocols['complex-explicit'] = alchemical_states

    # Create phases.
    phase = 'complex-explicit'
    systems = { phase : system }
    positions = { phase : positions }
    phases = [phase]
    atom_indices = { 'complex-explicit' : { 'ligand' : [0] } } # first particle is ligand

    # Create new simulation.
    yank = Yank(store_dir)
    yank.create(phases, systems, positions, atom_indices, thermodynamic_state, options=options, protocols=protocols)

    # Run the simulation.
    yank.run()

    # Analyze the data.
    results = yank.analyze()
    standard_state_correction = results[phase]['standard_state_correction']
    Delta_f = results[phase]['Delta_f_ij'][0,-1] + standard_state_correction
    dDelta_f = results[phase]['dDelta_f_ij'][0,-1]
    nsigma = abs(binding_free_energy/kT - Delta_f) / dDelta_f

    # Check results against analytical results.
    # TODO: Incorporate standard state correction
    output = ""
    output += "Analytical binding free energy                               : %10.5f +- %10.5f kT\n" % (binding_free_energy / kT, 0)
    output += "Computed binding free energy (with standard state correction): %10.5f +- %10.5f kT (nsigma = %3.1f)\n" % (Delta_f, dDelta_f, nsigma)
    output += "Standard state correction alone                              : %10.5f           kT\n" % (standard_state_correction)
    print output

    if (nsigma > NSIGMA_MAX):
        output += "\n"
        output += "Computed binding free energy differs from true binding free energy.\n"
        raise Exception(output)
<commit_msg>Minor change to standard-state correction direction<commit_after>#!/usr/local/bin/env python

"""
Test YANK using simple models.

DESCRIPTION

This test suite generates a number of simple models to test the 'Yank' facility.

COPYRIGHT

Written by <NAME> <<EMAIL>> while at the University of California Berkeley.

LICENSE

This code is licensed under the latest available version of the GNU General Public License.

"""

#=============================================================================================
# GLOBAL IMPORTS
#=============================================================================================

import os
import sys
import math
import copy
import time
import datetime
from functools import partial

import numpy

from simtk import openmm
from openmmtools import testsystems

from yank import Yank
from yank.repex import ThermodynamicState

import logging
logger = logging.getLogger(__name__)

#=============================================================================================
# MODULE CONSTANTS
#=============================================================================================

from simtk import unit
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA # Boltzmann constant

#=============================================================================================
# MAIN AND TESTS
#=============================================================================================

def test_LennardJonesPair():
    """
    Compute binding free energy of two Lennard-Jones particles and compare to numerical result.

    """

    NSIGMA_MAX = 6.0 # number of standard errors tolerated for success

    # Create Lennard-Jones pair.
    test = testsystems.LennardJonesPair()
    system, positions = test.system, test.positions
    thermodynamic_state = ThermodynamicState(temperature=300.0*unit.kelvin)
    binding_free_energy = test.get_binding_free_energy(thermodynamic_state)
    kT = kB * thermodynamic_state.temperature

    # Create temporary directory for testing.
    import tempfile
    store_dir = tempfile.mkdtemp()
    
    # Initialize YANK object.
    options = dict()
    options['restraint_type'] = 'flat-bottom'
    options['number_of_iterations'] = 50
    options['platform'] = openmm.Platform.getPlatformByName("Reference") # use Reference platform for speed
    options['mc_rotation'] = False
    options['mc_displacement'] = True
    options['mc_displacement_sigma'] = 1.0 * unit.nanometer
    options['verbose'] = False
    options['timestep'] = 5 * unit.femtoseconds
    options['nsteps_per_iteration'] = 20

    # Override box vectors.
    box_edge = 30.0 * unit.angstrom
    a = unit.Quantity((box_edge,        0 * unit.angstrom, 0 * unit.angstrom))
    b = unit.Quantity((0 * unit.angstrom, box_edge,        0 * unit.angstrom))
    c = unit.Quantity((0 * unit.angstrom, 0 * unit.angstrom, box_edge))
    system.setDefaultPeriodicBoxVectors(a, b, c)

    # Alchemical protocol.
    from yank.alchemy import AlchemicalState
    alchemical_states = list()
    lambda_values = [0.0, 0.25, 0.50, 0.75, 1.0]
    for lambda_value in lambda_values:
        alchemical_state = AlchemicalState()
        alchemical_state['lambda_electrostatics'] = lambda_value
        alchemical_state['lambda_sterics'] = lambda_value
        alchemical_states.append(alchemical_state)
    protocols = dict()
    protocols['complex-explicit'] = alchemical_states

    # Create phases.
    phase = 'complex-explicit'
    systems = { phase : system }
    positions = { phase : positions }
    phases = [phase]
    atom_indices = { 'complex-explicit' : { 'ligand' : [0] } } # first particle is ligand

    # Create new simulation.
    yank = Yank(store_dir)
    yank.create(phases, systems, positions, atom_indices, thermodynamic_state, options=options, protocols=protocols)

    # Run the simulation.
    yank.run()

    # Analyze the data.
    results = yank.analyze()
    standard_state_correction = results[phase]['standard_state_correction']
    Delta_f = results[phase]['Delta_f_ij'][0,-1] - standard_state_correction
    dDelta_f = results[phase]['dDelta_f_ij'][0,-1]
    nsigma = abs(binding_free_energy/kT - Delta_f) / dDelta_f

    # Check results against analytical results.
    # TODO: Incorporate standard state correction
    output = ""
    output += "Analytical binding free energy                               : %10.5f +- %10.5f kT\n" % (binding_free_energy / kT, 0)
    output += "Computed binding free energy (with standard state correction): %10.5f +- %10.5f kT (nsigma = %3.1f)\n" % (Delta_f, dDelta_f, nsigma)
    output += "Standard state correction alone                              : %10.5f           kT\n" % (standard_state_correction)
    print output

    if (nsigma > NSIGMA_MAX):
        output += "\n"
        output += "Computed binding free energy differs from true binding free energy.\n"
        raise Exception(output)

{'fraction_non_alphanumeric': 0.19953051643192488, 'fraction_numerical': 0.011839150847111655, 'mean_word_length': 3.994393476044852, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7355125', 'n_tokens_mistral': 2885, 'n_tokens_neox': 2662, 'n_words': 721}