Skip to content

Commit

Permalink
Merge pull request #172 from kyleabeauchamp/del
Browse files Browse the repository at this point in the history 
Deleted several more unused functions.
  • Loading branch information
@kyleabeauchamp
kyleabeauchamp committed Jan 12, 2015
2 parents 87d52e5 + 7748766 commit 1c0268e
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 385 deletions.
357 changes: 0 additions & 357 deletions Yank/analyze.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,292 +39,6 @@
#=============================================================================================
# SUBROUTINES
#=============================================================================================
def write_file(filename, contents):
"""Write the specified contents to a file.
ARGUMENTS
filename (string) - the file to be written
contents (string) - the contents of the file to be written
"""

outfile = open(filename, 'w')

if type(contents) == list:
for line in contents:
outfile.write(line)
elif type(contents) == str:
outfile.write(contents)
else:
raise "Type for 'contents' not supported: " + repr(type(contents))

outfile.close()

return

def read_file(filename):
"""Read contents of the specified file.
ARGUMENTS
filename (string) - the name of the file to be read
RETURNS
lines (list of strings) - the contents of the file, split by line
"""

infile = open(filename, 'r')
lines = infile.readlines()
infile.close()

return lines

def initialize_netcdf(netcdf_file, title, natoms, is_periodic = False, has_velocities = False):
"""Initialize the given NetCDF file according to the AMBER NetCDF Convention Version 1.0, Revision B.
ARGUMENTS
netcdf_file (NetCDFFile object) - the file to initialize global attributes, dimensions, and variables for
title (string) - the title for the netCDF file
natoms (integer) - the number of atoms in the trajectories to be written
is_periodic (boolean) - if True, box coordinates will also be stored
has_velocities (boolean) - if True, the velocity trajectory variables will also be created
NOTES
The AMBER NetCDF convention is defined here:
http://amber.scripps.edu/netcdf/nctraj.html
"""

# Create dimensions.
netcdf_file.createDimension('frame', 0) # unlimited number of frames in trajectory
netcdf_file.createDimension('spatial', 3) # number of spatial coordinates
netcdf_file.createDimension('atom', natoms) # number of atoms in the trajectory
netcdf_file.createDimension('label', 5) # label lengths for cell dimensions
netcdf_file.createDimension('cell_spatial', 3) # naming conventions for cell spatial dimensions
netcdf_file.createDimension('cell_angular', 3) # naming conventions for cell angular dimensions

# Set attributes.
setattr(netcdf_file, 'title', title)
setattr(netcdf_file, 'application', 'AMBER')
setattr(netcdf_file, 'program', 'sander')
setattr(netcdf_file, 'programVersion', '8')
setattr(netcdf_file, 'Conventions', 'AMBER')
setattr(netcdf_file, 'ConventionVersion', '1.0')

# Define variables to store unit cell data, if specified.
if is_periodic:
cell_spatial = netcdf_file.createVariable('cell_spatial', 'c', ('cell_spatial',))
cell_angular = netcdf_file.createVariable('cell_angular', 'c', ('cell_spatial', 'label'))
cell_lengths = netcdf_file.createVariable('cell_lengths', 'd', ('frame', 'cell_spatial'))
setattr(cell_lengths, 'units', 'angstrom')
cell_angles = netcdf_file.createVariable('cell_angles', 'd', ('frame', 'cell_angular'))
setattr(cell_angles, 'units', 'degree')

netcdf_file.variables['cell_spatial'][0] = 'x'
netcdf_file.variables['cell_spatial'][1] = 'y'
netcdf_file.variables['cell_spatial'][2] = 'z'

netcdf_file.variables['cell_angular'][0] = 'alpha'
netcdf_file.variables['cell_angular'][1] = 'beta '
netcdf_file.variables['cell_angular'][2] = 'gamma'

# Define variables to store velocity data, if specified.
if has_velocities:
velocities = netcdf_file.createVariable('velocities', 'd', ('frame', 'atom', 'spatial'))
setattr(velocities, 'units', 'angstrom/picosecond')
setattr(velocities, 'scale_factor', 20.455)

# Define coordinates and snapshot times.
frame_times = netcdf_file.createVariable('time', 'f', ('frame',))
setattr(frame_times, 'units', 'picosecond')
frame_coordinates = netcdf_file.createVariable('coordinates', 'f', ('frame', 'atom', 'spatial'))
setattr(frame_coordinates, 'units', 'angstrom')

# Define optional data not specified in the AMBER NetCDF Convention that we will make use of.
frame_energies = netcdf_file.createVariable('total_energy', 'f', ('frame',))
setattr(frame_energies, 'units', 'kilocalorie/mole')
frame_energies = netcdf_file.createVariable('potential_energy', 'f', ('frame',))
setattr(frame_energies, 'units', 'kilocalorie/mole')

return

def write_netcdf_frame(netcdf_file, frame_index, time = None, coordinates = None, cell_lengths = None, cell_angles = None, total_energy = None, potential_energy = None):
"""Write a NetCDF frame.
ARGUMENTS
netcdf_file (NetCDFFile) - the file to write a frame to
frame_index (integer) - the frame to be written
OPTIONAL ARGUMENTS
time (float) - time of frame (in picoseconds)
coordinates (natom x nspatial NumPy array) - atomic coordinates (in Angstroms)
cell_lengths (nspatial NumPy array) - cell lengths (Angstroms)
cell_angles (nspatial NumPy array) - cell angles (degrees)
total_energy (float) - total energy (kcal/mol)
potential_energy (float) - potential energy (kcal/mol)
"""
if time != None: netcdf_file.variables['time'][frame_index] = time
if coordinates != None: netcdf_file.variables['coordinates'][frame_index,:,:] = coordinates
if cell_lengths != None: netcdf_file.variables['cell_lengths'][frame_index,:] = cell_lengths
if cell_angles != None: netcdf_file.variables['cell_angles'][frame_index,:] = cell_angles
if total_energy != None: netcdf_file.variables['total_energy'][frame_index] = total_energy
if potential_energy != None: netcdf_file.variables['total_energy'][frame_index] = potential_energy

return



def write_netcdf_replica_trajectories(directory, prefix, title, ncfile):
"""Write out replica trajectories in AMBER NetCDF format.
ARGUMENTS
directory (string) - the directory to write files to
prefix (string) - prefix for replica trajectory files
title (string) - the title to give each NetCDF file
ncfile (NetCDF) - NetCDF file object for input file
"""
# Get current dimensions.
niterations = ncfile.variables['positions'].shape[0]
nstates = ncfile.variables['positions'].shape[1]
natoms = ncfile.variables['positions'].shape[2]

# Write out each replica to a separate file.
for replica in range(nstates):
# Create a new replica file.
output_filename = os.path.join(directory, '%s-%03d.nc' % (prefix, replica))
ncoutfile = netcdf.Dataset(output_filename, 'w')
initialize_netcdf(ncoutfile, title + " (replica %d)" % replica, natoms)
for iteration in range(niterations):
coordinates = np.array(ncfile.variables['positions'][iteration,replica,:,:])
coordinates *= 10.0 # convert nm to angstroms
write_netcdf_frame(ncoutfile, iteration, time = 1.0 * iteration, coordinates = coordinates)
ncoutfile.close()

return


def write_pdb_replica_trajectories(basepdb, directory, prefix, title, ncfile, trajectory_by_state=True):
"""Write out replica trajectories as multi-model PDB files.
ARGUMENTS
basepdb (string) - name of PDB file to read atom names and residue information from
directory (string) - the directory to write files to
prefix (string) - prefix for replica trajectory files
title (string) - the title to give each PDB file
ncfile (NetCDF) - NetCDF file object for input file
"""
niterations = ncfile.variables['positions'].shape[0]
nstates = ncfile.variables['positions'].shape[1]
natoms = ncfile.variables['positions'].shape[2]

atom_list=read_pdb(basepdb)
if (len(atom_list) != natoms):
logger.info ("Number of atoms in trajectory (%d) differs from number of atoms in reference PDB (%d)." % (natoms, len(atom_list)))
raise Exception

if trajectory_by_state:
for state_index in range(0,nstates):
logger.info("Working on state %d / %d" % (state_index,nstates))
file_name= "%s-%03d.pdb" % (prefix,state_index)
full_filename=directory+'/'+file_name
outfile = open(full_filename, 'w')
for iteration in range(niterations):
state_indices = ncfile.variables['states'][iteration,:]
replica_index = list(state_indices).index(state_index)
outfile.write('MODEL %4d\n' % (iteration+1))
write_pdb(atom_list,outfile,iteration,replica_index,title,ncfile,trajectory_by_state=True)
outfile.write('ENDMDL\n')

outfile.close()

else:
for replica_index in range(nstates):
logger.info( "Working on replica %d / %d" % (replica_index,nstates))
file_name="R-%s-%03d.pdb" % (prefix,replica_index)
full_filename=directory+'/'+file_name
outfile = open(full_filename, 'w')
for iteration in range(niterations):
outfile.write('MODEL %4d\n' % (iteration+1))
write_pdb(atom_list,outfile,iteration,replica_index,title,ncfile,trajectory_by_state=False)
outfile.write('ENDMDL\n')

outfile.close()

return

def read_pdb(filename):
"""
Read the contents of a PDB file.
ARGUMENTS
filename (string) - name of the file to be read
RETURNS
atoms (list of dict) - atoms[index] is a dict of fields for the ATOM residue
"""

# Read the PDB file into memory.
pdbfile = open(filename, 'r')

# Extract the ATOM entries.
# Format described here: http://bmerc-www.bu.edu/needle-doc/latest/atom-format.html
atoms = list()
for line in pdbfile:
if line[0:6] == "ATOM ":
# Parse line into fields.
atom = dict()
atom["serial"] = line[6:11]
atom["atom"] = line[12:16]
atom["altLoc"] = line[16:17]
atom["resName"] = line[17:20]
atom["chainID"] = line[21:22]
atom["Seqno"] = line[22:26]
atom["iCode"] = line[26:27]
atom["x"] = line[30:38]
atom["y"] = line[38:46]
atom["z"] = line[46:54]
atom["occupancy"] = line[54:60]
atom["tempFactor"] = line[60:66]
atoms.append(atom)

# Close PDB file.
pdbfile.close()

# Return dictionary of present residues.
return atoms

def write_pdb(atoms, filename, iteration, replica, title, ncfile,trajectory_by_state=True):
"""Write out replica trajectories as multi-model PDB files.
ARGUMENTS
atoms (list of dict) - parsed PDB file ATOM entries from read_pdb() - WILL BE CHANGED
filename (string) - name of PDB file to be written
title (string) - the title to give each PDB file
ncfile (NetCDF) - NetCDF file object for input file
"""

# Extract coordinates to be written.
coordinates = np.array(ncfile.variables['positions'][iteration,replica,:,:])
coordinates *= 10.0 # convert nm to angstroms

# Create file.
#outfile = open(filename, 'w')

# Write ATOM records.
for (index, atom) in enumerate(atoms):
atom["x"] = "%8.3f" % coordinates[index,0]
atom["y"] = "%8.3f" % coordinates[index,1]
atom["z"] = "%8.3f" % coordinates[index,2]
filename.write('ATOM %(serial)5s %(atom)4s%(altLoc)c%(resName)3s %(chainID)c%(Seqno)5s %(x)8s%(y)8s%(z)8s\n' % atom)

# Close file.
#outfile.close()

return


def show_mixing_statistics(ncfile, cutoff=0.05, nequil=0):
Expand Down Expand Up @@ -384,77 +98,6 @@ def show_mixing_statistics(ncfile, cutoff=0.05, nequil=0):

return

def analyze_acceptance_probabilities(ncfile, cutoff = 0.4):
"""Analyze acceptance probabilities.
ARGUMENTS
ncfile (NetCDF) - NetCDF file to be analyzed.
OPTIONAL ARGUMENTS
cutoff (float) - cutoff for showing acceptance probabilities as blank (default: 0.4)
"""

# Get current dimensions.
niterations = ncfile.variables['mixing'].shape[0]
nstates = ncfile.variables['mixing'].shape[1]

# Compute mean.
mixing = ncfile.variables['mixing'][:,:,:]
Pij = np.mean(mixing, 0)

# Write title.
print "Average state-to-state acceptance probabilities"
print "(Probabilities less than %(cutoff)f shown as blank.)" % vars()
print ""

# Write header.
print "%4s" % "",
for j in range(nstates):
print "%6d" % j,
print ""

# Write rows.
for i in range(nstates):
print "%4d" % i,
for j in range(nstates):
if Pij[i,j] > cutoff:
print "%6.3f" % Pij[i,j],
else:
print "%6s" % "",

print ""

return


def check_positions(ncfile):
"""Make sure no positions have gone 'nan'.
ARGUMENTS
ncfile (NetCDF) - NetCDF file object for input file
"""

# Get current dimensions.
niterations = ncfile.variables['positions'].shape[0]
nstates = ncfile.variables['positions'].shape[1]
natoms = ncfile.variables['positions'].shape[2]

# Compute torsion angles for each replica
for iteration in range(niterations):
for replica in range(nstates):
# Extract positions
positions = np.array(ncfile.variables['positions'][iteration,replica,:,:])
# Check for nan
if np.any(np.isnan(positions)):
# Nan found -- raise error
logger.info("Iteration %d, state %d - nan found in positions." % (iteration, replica))
# Report coordinates
for atom_index in range(natoms):
logger.info("%16.3f %16.3f %16.3f" % (positions[atom_index,0], positions[atom_index,1], positions[atom_index,2]))
if np.any(np.isnan(positions[atom_index,:])):
raise "nan detected in positions"

return

def estimate_free_energies(ncfile, ndiscard = 0, nuse = None):
"""Estimate free energies of all alchemical states.
Expand Down
Loading

0 comments on commit 1c0268e

Please sign in to comment.