KaneAI - World’s First E2E Software Testing Agent.Book a Demo Join TestMu Conference

How to use molecule_data method in molecule

Best Python code snippet using molecule_python

lines2xsec.py

Source:lines2xsec.py

1#!/usr/bin/env python2#3# Copyright (C) 2017 - Massachusetts Institute of Technology (MIT)4#5# This program is free software: you can redistribute it and/or modify6# it under the terms of the GNU General Public License as published by7# the Free Software Foundation, either version 3 of the License, or8# (at your option) any later version.9#10# This program is distributed in the hope that it will be useful,11# but WITHOUT ANY WARRANTY; without even the implied warranty of12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the13# GNU General Public License for more details.14#15# You should have received a copy of the GNU General Public License16# along with this program.  If not, see <http://www.gnu.org/licenses/>.17"""18A Re-worked version of the old line_list_to_cross_section_calculation used in 0.4-0.6 19to generated molecular cross sections. 20This is still a bandage patch and does not solve some of the core problems with the fix21Will need a complete new look to check for accuracy in version 0.822"""23import os24import hapi as hp25import numpy as np26from bisect import bisect27from numpy import complex128,int64,float64,exp28SI = False29# default values for intensity threshold30DefaultIntensityThreshold = 0. # cm*molec31# default value for omega wing32DefaultOmegaWing = None33# default value for omega wing in halfwidths (from center)34DefaultOmegaWingHW = 50. # cm-1    HOTW default35# Defining precision36__ComplexType__ = complex12837__IntegerType__ = int6438__FloatType__ = float6439if SI:40    cBolts = 1.3806503e-2341    cc = 2.99792458e842else: #CGS43    cBolts = 1.380648813E-16 # erg/K, CGS44    cc = 2.99792458e10 # cm/s, CGS45# No Bias46cMassMol = 1.66053873e-2747cZero = 0.48# Hapi dependencies49# this is too much work to isolate and handle50PROFILE_VOIGT = hp.PROFILE_VOIGT   51PYTIPS=hp.PYTIPS52# Volume concentration of all gas molecules at the pressure p and temperature T53def volumeConcentration(p,T):54    if SI:55        return (p/9.869233e-6)/(cBolts*T) # SI56    else:57        return (p/9.869233e-7)/(cBolts*T) # CGS58# temperature dependence for intencities (HITRAN)59def EnvironmentDependency_Intensity(LineIntensityRef,T,Tref,SigmaT,SigmaTref,60                                    LowerStateEnergy,LineCenter):61    const = __FloatType__(1.4388028496642257)  # check this unit62    ch = exp(-const*LowerStateEnergy/T)*(1-exp(-const*LineCenter/T))63    zn = exp(-const*LowerStateEnergy/Tref)*(1-exp(-const*LineCenter/Tref))64    LineIntensity = LineIntensityRef*SigmaTref/SigmaT*ch/zn65    return LineIntensity66def EnvironmentDependency_Gamma0(Gamma0_ref,T,Tref,p,pref,TempRatioPower):67    return Gamma0_ref*p/pref*(Tref/T)**TempRatioPower68def read_data(path,molecule,numin,numax,imin=-1,ctop=-1,direct=False,multi=False):69    """70    numin: minimum wavenumber71    numax: maximum wavenumber72    imax:  minimum intensity73    ctop:  Top intensity  74    75    Reading data from locally stored line list files 76    Can be constrained in wavenumber77    Returns a dictionary containing data from the line list file78    """79    if multi:80        81        M, I, wavenumber, intensity,gamma_air,gamma_self, elower, n_air,delta_air = [],[],[],[],[],[],[],[],[]82            83        for filename in path:84            85            front,back = filename.split("-")86            87            filemin = float(front.split("_")[-1])88            filemax = float(back.split("_")[0])89            90            if filemin > numax or filemax < numin:91                continue92            93            94            #print "Reading Data From %s"%(filename)95            96            97            f = open(filename).read().split("\n")98            99            100            for i in range(len(f)-1):101                102                wn = float(f[i][3:15])103                104                if wn < numin or wn >= numax:105                    continue106            107                M.append(         int(f[i][0:2]))108                I.append(         int(f[i][2]))109                wavenumber.append(float(f[i][3:15]))110                intensity.append( float(f[i][16:26]))111                gamma_air.append( float(f[i][35:40]))112                gamma_self.append(float(f[i][41:45]))113                elower.append(    float(f[i][46:56]))114                n_air.append(     float(f[i][56:59]))115                delta_air.append( float(f[i][60:68]))116                117        data = {"M":M, "I":I, "wavenumber":wavenumber, "intensity":intensity, 118                  "gamma_air":gamma_air, "gamma_self":gamma_self, "elower":elower,119                  "n_air":n_air, "delta_air":delta_air}120        121        return data122    123    else:124        125        if direct:126            print("Reading Data From %s"%(path))127            f = open(path).read().split("\n")128            129        else:130            #print "Reading Data From %s/%s.data"%(path,molecule)131            f = open("%s/%s.data"%(path,molecule)).read().split("\n")132        133        M, I, wavenumber, intensity,gamma_air,gamma_self, elower, n_air,delta_air = [],[],[],[],[],[],[],[],[]134        135        appending = False136        for i in range(len(f)-1):137            138            wn = float(f[i][3:15])139            140            if wn >= numin:141                appending = True142            if wn >= numax:143                break144            145            M.append(         int(f[i][0:2]))146            I.append(         int(f[i][2]))147            wavenumber.append(float(f[i][3:15]))148            intensity.append( float(f[i][16:26]))149            gamma_air.append( float(f[i][35:40]))150            gamma_self.append(float(f[i][41:45]))151            elower.append(    float(f[i][46:56]))152            n_air.append(     float(f[i][56:59]))153            delta_air.append( float(f[i][60:68]))154            155        data = {"M":M, "I":I, "wavenumber":wavenumber, "intensity":intensity, 156                  "gamma_air":gamma_air, "gamma_self":gamma_self, "elower":elower,157                  "n_air":n_air, "delta_air":delta_air}158        159        return data160def absorption_Voigt_calculation(molecule_data, component, gamma_name, P, T, numin, numax, step=0.1, cross_section=True, OmegaWing=DefaultOmegaWing,161                 IntensityThreshold=DefaultIntensityThreshold, OmegaWingHW=DefaultOmegaWingHW, partitionFunction=PYTIPS):162    """163    Considering both Doppler Broadening and Pressure Broadening164    165    calculating cross section for each molecule166    167    P: Pressure in Pa, although need to look into P. 1atm should be 101300. For simplicity we're using 100000... for now.168    169    """170    171    # Reference Conditions172    Tref = 296.173    Pref = 1.174    P = P/100000.0175    176    177    # Determining the length of the wavenumber array178    nu = molecule_data["wavenumber"]179    180    181    Omega_min = float(numin)182    Omega_max = float(numax)183    OmegaStep = step   184    OmegaCount = (Omega_max-Omega_min)/OmegaStep+1185    Omegas = np.linspace(Omega_min,Omega_max,OmegaCount)[:-1]186    Xsect  = np.zeros(len(Omegas)) 187    188    189    ABUNDANCES = {}190    NATURAL_ABUNDANCES = {}191    192    M = molecule_data["M"][0]193    I = molecule_data["I"][0]194    if len(component) >=3:195        N = component[2]196    else:197        N = hp.ISO[(M,I)][hp.ISO_INDEX["abundance"]]198    199    ABUNDANCES[(M,I)] = N200    NATURAL_ABUNDANCES[(M,I)] = hp.ISO[(M,I)][hp.ISO_INDEX["abundance"]]    201    202    factor = hp.volumeConcentration(P,T)203    EnvDependences = lambda ENV, LINE:{}204    Env = {}205    Env["T"] = float(T)206    Env["P"] = float(P)207    Env["Tref"] = Tref208    Env["Pref"] = Pref    209    210        211    nline = len(molecule_data["wavenumber"])212    MolNumberDB   = int(molecule_data["M"][0])213    IsoNumberDB   = int(molecule_data["I"][0])214    215    if SI:216        m = hp.ISO[(MolNumberDB,IsoNumberDB)][hp.ISO_INDEX['mass']]*cMassMol* cMassMol 217    else:218        m = hp.ISO[(MolNumberDB,IsoNumberDB)][hp.ISO_INDEX['mass']]*cMassMol* 1000.219    220    parnames = molecule_data.keys()221    222    # alot of calculations can be omitted here given we don't need to calculate the entire thing223    for RowID in range(nline):224        225        Line = {parname:molecule_data[parname][RowID] for parname in parnames}226        CustomEnvDependences = EnvDependences(Env,Line)227        228        if (MolNumberDB,IsoNumberDB) not in ABUNDANCES: continue229        230        SigmaT = partitionFunction(MolNumberDB,IsoNumberDB,T)231        SigmaTref = partitionFunction(MolNumberDB,IsoNumberDB,Tref)   232         233        LineCenterDB = molecule_data['wavenumber'][RowID]234        235        236        LineIntensityDB = molecule_data['intensity'][RowID]237        LowerStateEnergyDB = molecule_data['elower'][RowID]238        if 'intensity' in CustomEnvDependences:239            LineIntensity = CustomEnvDependences['intensity']240        else:241            LineIntensity = EnvironmentDependency_Intensity(LineIntensityDB,T,Tref,SigmaT,SigmaTref,242                                                            LowerStateEnergyDB,LineCenterDB)243        244        #   FILTER by LineIntensity: compare it with IntencityThreshold245        if LineIntensity < IntensityThreshold: continue246        # Doppler Broadening247        GammaD = np.sqrt(2*cBolts*T*np.log(2)/m/cc**2)*LineCenterDB248        249        # Pressure Broadening250        Gamma0 = 0.251        Shift0 = 0.252        # default is air, not considering multiple diluent atm253        abun = 1254        try:255            Gamma0DB =  molecule_data[gamma_name][RowID]256        except:257            Gamma0DB = 0.0258        259        260        if gamma_name == "gamma_air":261            n_name = "n_air"262            d_name = "delta_air"263            dp_name= "deltap_air"264        elif gamma_name == "gamma_self":265            n_name = "n_self"266            d_name = "delta_self"267            dp_name= "deltap_self"268        else:269            print("unknown gamma")270            return271        272        try:273            TempRatioPowerDB = molecule_data[n_name][RowID]274            if n_name == "n_self" and TempRatioPowerDB == 0.:275                TempRatioPowerDB = molecule_data["n_air"][RowID]276        except:277            TempRatioPowerDB = molecule_data["n_air"][RowID]278        279        Gamma0 += abun*CustomEnvDependences.get(gamma_name,280                  EnvironmentDependency_Gamma0(Gamma0DB,T,Tref,P,Pref,TempRatioPowerDB))281        try:282            Shift0DB = molecule_data[d_name][RowID]283        except:284            Shift0DB = 0.0285        286        try:287            deltap = molecule_data[dp_name][RowID]288        except:289            deltap = 0.0290        291        # pressure broadening for self is 0?292        Shift0 += abun*CustomEnvDependences.get(d_name, # default ->293                  ((Shift0DB + deltap*(T-Tref))*P/Pref))294        #Shift0 = 0295    296        OmegaWingF = max(OmegaWing,OmegaWingHW*Gamma0,OmegaWingHW*GammaD)297        BoundIndexLower = bisect(Omegas,LineCenterDB-OmegaWingF)298        BoundIndexUpper = bisect(Omegas,LineCenterDB+OmegaWingF)299        300        lineshape_vals = PROFILE_VOIGT(LineCenterDB+Shift0,GammaD,Gamma0,Omegas[BoundIndexLower:BoundIndexUpper])[0]301        302        303        # absorption coefficient calculation304        305        306        if cross_section:307            Xsect[BoundIndexLower:BoundIndexUpper] += ABUNDANCES[(MolNumberDB,IsoNumberDB)]/NATURAL_ABUNDANCES[(MolNumberDB,IsoNumberDB)] * \308                                                      LineIntensity * lineshape_vals309        else:310            Xsect[BoundIndexLower:BoundIndexUpper] += factor / NATURAL_ABUNDANCES[(MolNumberDB,IsoNumberDB)] * \311                                                      ABUNDANCES[(MolNumberDB,IsoNumberDB)] * \312                                                      LineIntensity * lineshape_vals313        ...

analyzeall.py

Source:analyzeall.py

1#!/opt/local/bin/python2.52#=============================================================================================3# Analyze a set of datafiles produced by YANK.4#=============================================================================================5#=============================================================================================6# REQUIREMENTS7#8# The netcdf4-python module is now used to provide netCDF v4 support:9# http://code.google.com/p/netcdf4-python/10#11# This requires NetCDF with version 4 and multithreading support, as well as HDF5.12#=============================================================================================13#=============================================================================================14# TODO15#=============================================================================================16#=============================================================================================17# CHAGELOG18#=============================================================================================19#=============================================================================================20# VERSION CONTROL INFORMATION21#=============================================================================================22#=============================================================================================23# IMPORTS24#=============================================================================================25import numpy26from numpy import *27#from Scientific.IO import NetCDF # scientific python28#import scipy.io.netcdf as netcdf29import netCDF4 as netcdf # netcdf4-python30import os31import sys32import os.path33import math34import gzip35from pymbar import MBAR # multistate Bennett acceptance ratio36import timeseries # for statistical inefficiency analysis37import simtk.unit as units38from yank.analysis import *39#=============================================================================================40# MAIN41#=============================================================================================42#data_directory = '/Users/yank/data-from-lincoln/T4-lysozyme-L99A/amber-gbsa/amber-gbsa/' # directory containing datafiles43#data_directory = '/Users/yank/data-from-lincoln/FKBP/amber-gbsa/' # directory containing datafiles44#data_directory = '/Users/yank/data-from-lincoln/FKBP/amber-gbvi/' # directory containing datafiles45#data_directory = '/uf/ac/jchodera/code/yank/test-systems/T4-lysozyme-L99A/amber-gbsa/amber-gbsa/' # directory containing datafiles46#data_directory = '/scratch/users/jchodera/yank/test-systems/T4-lysozyme-L99A/amber-gbsa/amber-gbsa/' # directory containing datafiles47data_directory = 'examples/p-xylene' # directory containing datafiles48# Store molecule data.49molecule_data = dict()50# Generate list of files in this directory.51import commands52molecules = commands.getoutput('ls -1 %s' % data_directory).split()53for molecule in molecules:54    source_directory = os.path.join(data_directory, molecule)55    print source_directory56    # Storage for different phases.57    data = dict()58    phases = ['vacuum', 'solvent', 'complex']    59    # Process each netcdf file.60    for phase in phases:61        # Construct full path to NetCDF file.62        fullpath = os.path.join(source_directory, phase + '.nc')63        # Skip if the file doesn't exist.64        if (not os.path.exists(fullpath)): continue65        # Open NetCDF file for reading.66        print "Opening NetCDF trajectory file '%(fullpath)s' for reading..." % vars()67        ncfile = netcdf.Dataset(fullpath, 'r')68        # DEBUG69        print "dimensions:"70        for dimension_name in ncfile.dimensions.keys():71            print "%16s %8d" % (dimension_name, len(ncfile.dimensions[dimension_name]))72    73        # Read dimensions.74        niterations = ncfile.variables['positions'].shape[0]75        nstates = ncfile.variables['positions'].shape[1]76        natoms = ncfile.variables['positions'].shape[2]77        print "Read %(niterations)d iterations, %(nstates)d states" % vars()78#        # Compute torsion trajectories79#        if phase in ['complex', 'receptor']:80#            print "Computing torsions..."81#            compute_torsion_trajectories(ncfile, os.path.join(source_directory, phase + ".val111"))82#        # Write out replica trajectories83#        print "Writing replica trajectories...\n"84#        title = 'Source %(source_directory)s phase %(phase)s' % vars()        85#        write_netcdf_replica_trajectories(source_directory, phase, title, ncfile)86        # Read reference PDB file.87        if phase in ['vacuum', 'solvent']:88            reference_pdb_filename = os.path.join(source_directory, "ligand.pdb")89        else:90            reference_pdb_filename = os.path.join(source_directory, "complex.pdb")91        atoms = read_pdb(reference_pdb_filename)92        # Write replica trajectories.93        #title = 'title'94        #write_pdb_replica_trajectories(reference_pdb_filename, source_directory, phase, title, ncfile, trajectory_by_state=False)95        96        # Check to make sure no self-energies go nan.97        check_energies(ncfile, atoms)98        # Check to make sure no positions are nan99        check_positions(ncfile)100        # Choose number of samples to discard to equilibration101        #nequil = 50102        #if phase == 'complex':103        #    nequil = 2000 # discard 2 ns of complex simulations104        u_n = extract_u_n(ncfile)105        [nequil, g_t, Neff_max] = detect_equilibration(u_n)106        print [nequil, Neff_max]107 108        # Examine acceptance probabilities.109        show_mixing_statistics(ncfile, cutoff=0.05, nequil=nequil)110        # Estimate free energies.111        (Deltaf_ij, dDeltaf_ij) = estimate_free_energies(ncfile, ndiscard = nequil)112    113        # Estimate average enthalpies114        (DeltaH_i, dDeltaH_i) = estimate_enthalpies(ncfile, ndiscard = nequil)115    116        # Accumulate free energy differences117        entry = dict()118        entry['DeltaF'] = Deltaf_ij[0,nstates-1] 119        entry['dDeltaF'] = dDeltaf_ij[0,nstates-1]120        entry['DeltaH'] = DeltaH_i[nstates-1] - DeltaH_i[0]121        entry['dDeltaH'] = numpy.sqrt(dDeltaH_i[0]**2 + dDeltaH_i[nstates-1]**2)122        data[phase] = entry123        # Get temperatures.124        ncvar = ncfile.groups['thermodynamic_states'].variables['temperatures']125        temperature = ncvar[0] * units.kelvin126        kT = kB * temperature127        # Close input NetCDF file.128        ncfile.close()129    # Skip if we have no data.130    if not ('vacuum' in data) or ('solvent' in data) or ('complex' in data): continue131    132    if (data.haskey('vacuum') and data.haskey('solvent')):133        # Compute hydration free energy (free energy of transfer from vacuum to water)134        DeltaF = data['vacuum']['DeltaF'] - data['solvent']['DeltaF']135        dDeltaF = numpy.sqrt(data['vacuum']['dDeltaF']**2 + data['solvent']['dDeltaF']**2)136        print "Hydration free energy: %.3f +- %.3f kT (%.3f +- %.3f kcal/mol)" % (DeltaF, dDeltaF, DeltaF * kT / units.kilocalories_per_mole, dDeltaF * kT / units.kilocalories_per_mole)137        # Compute enthalpy of transfer from vacuum to water138        DeltaH = data['vacuum']['DeltaH'] - data['solvent']['DeltaH']139        dDeltaH = numpy.sqrt(data['vacuum']['dDeltaH']**2 + data['solvent']['dDeltaH']**2)140        print "Enthalpy of hydration: %.3f +- %.3f kT (%.3f +- %.3f kcal/mol)" % (DeltaH, dDeltaH, DeltaH * kT / units.kilocalories_per_mole, dDeltaH * kT / units.kilocalories_per_mole)141    # Read standard state correction free energy.142    DeltaF_restraints = 0.0143    phase = 'complex'144    fullpath = os.path.join(source_directory, phase + '.nc')145    ncfile = netcdf.Dataset(fullpath, 'r')146    DeltaF_restraints = ncfile.groups['metadata'].variables['standard_state_correction'][0]147    ncfile.close()148    149    # Compute binding free energy (free energy of transfer from vacuum to water)150    DeltaF = data['solvent']['DeltaF'] - DeltaF_restraints - data['complex']['DeltaF']151    dDeltaF = numpy.sqrt(data['solvent']['dDeltaF']**2 + data['complex']['dDeltaF']**2)152    print ""153    print "Binding free energy : %16.3f +- %.3f kT (%16.3f +- %.3f kcal/mol)" % (DeltaF, dDeltaF, DeltaF * kT / units.kilocalories_per_mole, dDeltaF * kT / units.kilocalories_per_mole)154    print ""155    print "DeltaG vacuum       : %16.3f +- %.3f kT" % (data['vacuum']['DeltaF'], data['vacuum']['dDeltaF'])156    print "DeltaG solvent      : %16.3f +- %.3f kT" % (data['solvent']['DeltaF'], data['solvent']['dDeltaF'])157    print "DeltaG complex      : %16.3f +- %.3f kT" % (data['complex']['DeltaF'], data['complex']['dDeltaF'])158    print "DeltaG restraint    : %16.3f          kT" % DeltaF_restraints159    print ""160    # Compute binding enthalpy161    DeltaH = data['solvent']['DeltaH'] - DeltaF_restraints - data['complex']['DeltaH'] 162    dDeltaH = numpy.sqrt(data['solvent']['dDeltaH']**2 + data['complex']['dDeltaH']**2)163    print "Binding enthalpy    : %16.3f +- %.3f kT (%16.3f +- %.3f kcal/mol)" % (DeltaH, dDeltaH, DeltaH * kT / units.kilocalories_per_mole, dDeltaH * kT / units.kilocalories_per_mole)164    # Store molecule data.165    molecule_data[molecule] = data166# Extract sorted binding affinities.167sorted_molecules = ['1-methylpyrrole',168                    '1,2-dichlorobenzene',169                    '2-fluorobenzaldehyde',170                    '2,3-benzofuran',171                    'benzene',172                    'ethylbenzene',173                    'indene',174                    'indole',175                    'isobutylbenzene',176                    'n-butylbenzene',177                    'N-methylaniline',178                    'n-propylbenzene',179                    'o-xylene',180                    'p-xylene',181                    'phenol',182                    'toluene']183print ""184print "DeltaG"                                                                           185for molecule in sorted_molecules:186    try:187        DeltaF = molecule_data[molecule]['solvent']['DeltaF'] - molecule_data[molecule]['DeltaF_restraints'] - molecule_data[molecule]['complex']['DeltaF']188        dDeltaF = sqrt(molecule_data[molecule]['solvent']['dDeltaF']**2 + molecule_data[molecule]['complex']['dDeltaF']**2)189        print "%8.3f %8.3f %% %s" % (DeltaF, dDeltaF, molecule)190    except:191        print "%8.3f %8.3f %% %s" % (0.0, 0.0, molecule)        192        pass193print ""194print "DeltaH"                                                                           195for molecule in sorted_molecules:196    try:197        DeltaH = molecule_data[molecule]['solvent']['DeltaH'] - molecule_data[molecule]['complex']['DeltaH']198        dDeltaH = sqrt(molecule_data[molecule]['solvent']['dDeltaH']**2 + molecule_data[molecule]['complex']['dDeltaH']**2)199        print "%8.3f %8.3f %% %s" % (DeltaH, dDeltaH, molecule)200    except:201        print "%8.3f %8.3f %% %s" % (0.0, 0.0, molecule)                ...

rdkit_tutorial.py

Source:rdkit_tutorial.py

1import numpy as np2from rdkit import Chem3# Creating a simple molecule4toluene = Chem.MolFromSmiles('C1C=CC=CC1C')5hydrogen_cyanide = Chem.MolFromSmiles('C#N')6hydrogen_cyanide = Chem.MolFromSmiles('C#N')7# Build bond array8molecule = toluene9molecule_data = []10for i in range(toluene.GetNumAtoms()):11    molecule_data.append([])12    for j in range(toluene.GetNumAtoms()):13        if i == j or toluene.GetBondBetweenAtoms(i,j) == None:14            molecule_data[i].append(0)15        elif toluene.GetBondBetweenAtoms(i,j).GetBondType() == Chem.rdchem.BondType.SINGLE:16            molecule_data[i].append(1)17        elif toluene.GetBondBetweenAtoms(i,j).GetBondType() == Chem.rdchem.BondType.DOUBLE:18            molecule_data[i].append(2)19        elif toluene.GetBondBetweenAtoms(i, j).GetBondType() == Chem.rdchem.BondType.TRIPLE:20            molecule_data[i].append(3)21        else:22            molecule_data[i].append(0)...

Automation Testing Tutorials

Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.