Updated documentation, revised pycharmm engine

Author: Ryan

INSTALL

@@ -123,19 +123,23 @@ python -c "import alf"
 
 Installing without CUDA
 
-If you don't have access to CUDA GPUs, you can still run ALF, it will just be
-slower, and the tools aren't as streamlined for this. When setting up your
-environment don't include cuda. When you compile alf/wham, the lack of a cuda
-compiler should signal to cmake to compile the CPU code instead. This code has
-not been tested for a while, so please contact the developers if you have problems compiling. Next, you can probably skip compiling alf/dca unless you plan to
-estimate free energies with the Potts estimator. If you want to estimate free
-energies with the Potts estimator, you will need to modify CMakeLists.txt to
-omit plm. lm uses likelihood maximization which is better for most chemical
-spaces, and plm uses pseudolikelihood maximization which is typically only
-required for massive chemical spaces, e.g. more than a million species.
-
-You will need to modify your CHARMM input scripts to not use GPUs (see the
-scripts in alf/default_scripts). Removing the BLaDE or DOMDEC options, or
-replacing "domdec gpu only" with "domdec gpu off" in the scripts should work.
-You may wish to edit the scripts and the CHARMM calls in alf/runflat.py and
-alf/runprod.py to use multiple CPUs in parallel to improve efficiency.
+If you don't have access to CUDA GPUs, you can still run ALF, it will
+just be slower, and the tools aren't as streamlined for this. When
+setting up your environment don't include cuda. When you compile
+alf/wham, the lack of a cuda compiler should signal to cmake to compile
+the CPU code instead. This code has not been tested for a while, so
+please contact the developers if you have problems compiling. Next, you
+can probably skip compiling alf/dca unless you plan to estimate free
+energies with the Potts estimator. If you want to estimate free energies
+with the Potts estimator, you will need to modify CMakeLists.txt to omit
+plm. lm uses likelihood maximization which is better for most chemical
+spaces, and plm uses pseudolikelihood maximization which is typically
+only required for massive chemical spaces, e.g. more than a million
+species.
+
+You will need to modify your CHARMM input scripts to not use GPUs (see
+the scripts in alf/default_scripts). Removing the BLaDE or DOMDEC
+options, or replacing "domdec gpu only" with "domdec gpu off" in the
+scripts should work. You may wish to edit the scripts and the CHARMM
+calls in alf/runflat.py and alf/runprod.py to use multiple CPUs in
+parallel to improve efficiency.

README.md

@@ -36,11 +36,18 @@ Molecular dynamics in the alf software is performed by an external
 molecular dynamics engine which must be compiled and installed
 independently. The currently supported engines are:
  * charmm : the CHARMM software package utilizing the DOMDEC GPU ONLY
-   command
+   command for GPU acceleration
  * bladelib : the CHARMM software package utilizing the BLaDE library
+   for GPU acceleration. BLaDE is faster that DOMDEC, but has fewer
+   features
  * blade : the standalone BLaDE software package
+ * pycharmm : the python pyCHARMM package using CHARMM as a library to
+   call BLaDE
 These engines may be passed to alf routines to specify the engine,
 because some routines have slight differences based on the engine used.
+Use of other molecular dynamics engines in CHARMM is possible, including
+engines that do not use GPUs, see instructions at the end of INSTALL.
+Use of other engines will require more extensive modification of alf.
 
 Sampling will be optimal when the free energy landscape is flat, so
 flattening seeks to identify the biases that give flat free energy
@@ -168,11 +175,11 @@ higher order terms are assumed to be zero. This approximation is
 reasonable in many systems and significantly reduces sampling
 requirements, making it possible to estimate free energies for tens of
 thousands of sequences. To use the Potts model estimator, run
-alf.postprocess first, then see the examples in examples/engines with
-the _withPotts suffix. SubsetLM.sh uses likelihood maximization, best
-for systems with less than a million chemical end states, and
-SubsetPLM.sh uses pseudolikelihood maximization, best for systems with
-more than a million chemical end states. These use routines defined in
+alf.postprocess first, then see the examples in examples/engines/PottsLM
+and examples/engines/PottsPLM. PottsLM uses likelihood maximization,
+best for systems with less than a million chemical end states, and
+PottsPLM uses pseudolikelihood maximization, best for systems with more
+than a million chemical end states. These use routines defined in
 alf/dca.
 
 Coupling between sites: if you have multiple sites, it is possible these

alf/GetLambdas.py

@@ -1,6 +1,41 @@
 #! /usr/bin/env python
 
 def GetLambdas(alf_info,istep,ndupl=None,begres=None,endres=None):
+  """
+  Reads alchemical trajectories from binary format
+
+  This routine reads binary alchemical flattening trajectories from
+  run[i]/res/[name]_flat.lmd or binary alchemical production trajectories
+  from run[i][a]/res/[name]_prod[itt].lmd where [i] is the cycle number,
+  [a] is the duplicate letter, [name] is the system name, and [itt] is the
+  production chunk, and copies them into human readable trajectories in
+  analysis[i]/data/Lambda.[ia].[ir].dat where [ia] is the duplicate index
+  and [ir] is the replica index. This routine should be called from the
+  analysis[i] directory.
+
+  This routine can be called during flattening or production. Flattening
+  versus production is detected by the absence or presence, respectively
+  of the three optional parameters.
+
+  Parameters
+  ----------
+  alf_info : dict
+      Dictionary of variables alf needs to run
+  istep : int
+      The current cycle of alf being analyzed
+  ndupl : int, optional
+      The number of independent trials run in production. Leave empty to
+      signal this is flattening. (defaul is None)
+WORKING
+  Ff : int
+      The final cycle of alf to include in analysis (inclusive)
+  skipE : int, optional
+      In longer production runs the number of lambda samples may require
+      significant amounts of memory to store and analyze. Only alchemical
+      frames with index modulus skipE equal to skipE-1 are analyzed.
+      (default is 1 to analyze all frames) 
+  """
+
   import sys, os
   import numpy as np
   # from subprocess import call

alf/GetVolumes.py

@@ -27,7 +27,7 @@ def GetVolumes(alf_info,istep,ndupl=None,begres=None,endres=None):
   name=alf_info['name']
 
   if not 'q' in alf_info:
-    print("No charge list 'q' in alf_info")
+    print("No charge list 'q' in alf_info - no charge changing correction will be applied")
     return
   else:
     q=alf_info['q']

alf/SetVars.py

@@ -11,7 +11,6 @@ def SetVarsCharmm(alf_info,Step,minimize=False):
   nnodes=alf_info['nnodes']
   temp=alf_info['temp']
 
-
   fp=open('../variables'+str(Step)+'.inp','w')
   fp.write("* Variables from step %d of ALF\n*\n\n" % (Step,))
 
@@ -97,7 +96,7 @@ def SetVarsCharmm(alf_info,Step,minimize=False):
   for i in range(0,len(nsubs)):
     fp.write("set nsubs"+str(i+1)+" = "+str(nsubs[i])+"\n")
   fp.write("set temp = "+str(temp)+"\n")
-  fp.write("set minimize = "+str(int(minimize==True))+"\n")
+  fp.write("set minimizeflag = "+str(int(minimize==True))+"\n")
   fp.write("\n")
   fp.close()
 
@@ -114,7 +113,6 @@ def SetVarsBlade(alf_info,Step,minimize=False):
   nnodes=alf_info['nnodes']
   temp=alf_info['temp']
 
-
   fp=open('../variables'+str(Step)+'.inp','w')
 
   b_prev=np.loadtxt('b_prev.dat')
@@ -209,44 +207,80 @@ def SetVarsPycharmm(alf_info,Step,minimize=False):
   import yaml
   import copy
 
+  nblocks=alf_info['nblocks']
+  nsubs=alf_info['nsubs']
+  nreps=alf_info['nreps']
+  ncentral=alf_info['ncentral']
+  name=alf_info['name']
+  nnodes=alf_info['nnodes']
+  temp=alf_info['temp']
+
   fp=open('../variables'+str(Step)+'.py','w')
 
   fp.write("import yaml\n")
   fp.write("import numpy as np\n")
 
   bias={}
 
+  sub0=np.cumsum(nsubs)-nsubs
+
   b_prev=np.loadtxt('b_prev.dat')
   b=np.loadtxt('b.dat')
   b_sum=b_prev+b
   b_sum=np.reshape(b_sum,(1,-1))
   np.round(b_sum,decimals=2)
   np.savetxt('b_sum.dat',b_sum,fmt=' %7.2f')
 
-  bias['b']=b_sum.tolist()
+  for i in range(0,len(nsubs)):
+    for j in range(0,nsubs[i]):
+      key=f'lams{i+1}s{j+1}'
+      bias[key]=b_sum[0,sub0[i]+j].tolist()
 
   c_prev=np.loadtxt('c_prev.dat')
   c=np.loadtxt('c.dat')
   c_sum=c_prev+c
   np.round(c_sum,decimals=2)
   np.savetxt('c_sum.dat',c_sum,fmt=' %7.2f')
 
-  bias['c']=c_sum.tolist()
+  for si in range(0,len(nsubs)):
+    for sj in range(si,len(nsubs)):
+      for i in range(0,nsubs[si]):
+        j0=(i+1 if si==sj else 0)
+        for j in range(j0,nsubs[sj]):
+          key=f'cs{si+1}s{i+1}s{sj+1}s{j+1}'
+          bias[key]=-c_sum[sub0[si]+i,sub0[sj]+j].tolist()
 
   x_prev=np.loadtxt('x_prev.dat')
   x=np.loadtxt('x.dat')
   x_sum=x_prev+x
   np.round(x_sum,decimals=2)
   np.savetxt('x_sum.dat',x_sum,fmt=' %7.2f')
 
-  bias['x']=x_sum.tolist()
+  for si in range(0,len(nsubs)):
+    for sj in range(0,len(nsubs)):
+      for i in range(0,nsubs[si]):
+        for j in range(0,nsubs[sj]):
+          if sub0[si]+i!=sub0[sj]+j:
+            key=f'xs{si+1}s{i+1}s{sj+1}s{j+1}'
+            bias[key]=-x_sum[sub0[si]+i,sub0[sj]+j].tolist()
 
   s_prev=np.loadtxt('s_prev.dat')
   s=np.loadtxt('s.dat')
   s_sum=s_prev+s
   np.round(s_sum,decimals=2)
   np.savetxt('s_sum.dat',s_sum,fmt=' %7.2f')
 
+  for si in range(0,len(nsubs)):
+    for sj in range(0,len(nsubs)):
+      for i in range(0,nsubs[si]):
+        for j in range(0,nsubs[sj]):
+          if sub0[si]+i!=sub0[sj]+j:
+            key=f'ss{si+1}s{i+1}s{sj+1}s{j+1}'
+            bias[key]=-s_sum[sub0[si]+i,sub0[sj]+j].tolist()
+
+  bias['b']=b_sum.tolist()
+  bias['c']=c_sum.tolist()
+  bias['x']=x_sum.tolist()
   bias['s']=s_sum.tolist()
 
   fp.write("bias_string=\"\"\"\n")
@@ -258,86 +292,6 @@ def SetVarsPycharmm(alf_info,Step,minimize=False):
   fp.write("bias['x']=np.array(bias['x'])\n")
   fp.write("bias['s']=np.array(bias['s'])\n")
 
-  nsubs=alf_info['nsubs']
-
-  ibuff=0
-  lamss={}
-  for i in range(0,len(nsubs)):
-    lamss[i]={}
-    for j in range(0,nsubs[i]):
-      lamss[i][j]=b_sum[0,ibuff+j].tolist()
-    ibuff+=nsubs[i]
-
-  ibuff=0
-  cssss={}
-  for si in range(0,len(nsubs)):
-    jbuff=ibuff
-    cssss[si]={}
-    for sj in range(si,len(nsubs)):
-      cssss[si][sj]={}
-      for i in range(0,nsubs[si]):
-        ii=i+ibuff
-        cssss[si][sj][i]={}
-        j0=0
-        if si==sj:
-          j0=i+1
-        for j in range(j0,nsubs[sj]):
-          jj=j+jbuff
-          cssss[si][sj][i][j]=-c_sum[ii,jj].tolist()
-      jbuff+=nsubs[sj]
-    ibuff+=nsubs[si]
-
-  ibuff=0
-  xssss={}
-  for si in range(0,len(nsubs)):
-    xssss[si]={}
-    jbuff=0
-    for sj in range(0,len(nsubs)):
-      xssss[si][sj]={}
-      for i in range(0,nsubs[si]):
-        ii=i+ibuff
-        xssss[si][sj][i]={}
-        for j in range(0,nsubs[sj]):
-          jj=j+jbuff
-          if ii!=jj:
-            xssss[si][sj][i][j]=-x_sum[ii,jj].tolist()
-      jbuff+=nsubs[sj]
-    ibuff+=nsubs[si]
-
-  ibuff=0
-  sssss={}
-  for si in range(0,len(nsubs)):
-    sssss[si]={}
-    jbuff=0
-    for sj in range(0,len(nsubs)):
-      sssss[si][sj]={}
-      for i in range(0,nsubs[si]):
-        ii=i+ibuff
-        sssss[si][sj][i]={}
-        for j in range(0,nsubs[sj]):
-          jj=j+jbuff
-          if ii!=jj:
-            sssss[si][sj][i][j]=-s_sum[ii,jj].tolist()
-      jbuff+=nsubs[sj]
-    ibuff+=nsubs[si]
-
-  fp.write("lamss_string=\"\"\"\n")
-  yaml.dump(lamss,fp)
-  fp.write("\"\"\"\n")
-  fp.write("lamss=yaml.load(lamss_string,Loader=yaml.Loader)\n")
-  fp.write("cssss_string=\"\"\"\n")
-  yaml.dump(cssss,fp)
-  fp.write("\"\"\"\n")
-  fp.write("cssss=yaml.load(cssss_string,Loader=yaml.Loader)\n")
-  fp.write("xssss_string=\"\"\"\n")
-  yaml.dump(xssss,fp)
-  fp.write("\"\"\"\n")
-  fp.write("xssss=yaml.load(xssss_string,Loader=yaml.Loader)\n")
-  fp.write("sssss_string=\"\"\"\n")
-  yaml.dump(sssss,fp)
-  fp.write("\"\"\"\n")
-  fp.write("sssss=yaml.load(sssss_string,Loader=yaml.Loader)\n")
-
   alf_info_copy=copy.deepcopy(alf_info)
   alf_info_copy['nsubs']=alf_info_copy['nsubs'].tolist()
   alf_info_copy['nblocks']=alf_info_copy['nblocks'].tolist()
@@ -352,6 +306,17 @@ def SetVarsPycharmm(alf_info,Step,minimize=False):
   fp.write("if 'q' in alf_info:\n")
   fp.write("  alf_info['q']=np.array(alf_info['q'])\n")
 
+  fp.write("sysname='"+name+"'\n")
+  fp.write("nnodes="+str(nnodes)+"\n")
+  fp.write("nreps="+str(nreps)+"\n")
+  fp.write("ncentral="+str(ncentral)+"\n")
+  fp.write("nblocks="+str(nblocks)+"\n")
+  fp.write("nsites="+str(len(nsubs))+"\n")
+  fp.write("nsubs="+str(nsubs)+"\n")
+  for i in range(0,len(nsubs)):
+    fp.write("nsubs"+str(i+1)+"="+str(nsubs[i])+"\n")
+  fp.write("temp="+str(temp)+"\n")
+
   if minimize==True:
     fp.write("minimizeflag=True\n")
   else:

alf/default_scripts/blade_flat.inp

@@ -1,6 +1,14 @@
 ! BLaDE input file for lambda dynamics
 !
 
+! This script expects esteps (equilibration steps) and nsteps
+! (production steps) to be set by the calling process in
+! arguments.inp
+!
+! Other important variables such as temp (the system temperature), and
+! sysname (the system name in the prep directory) are set by
+! variablesflat.inp. minimizeflag is unused in blade
+
 verbose 0
 
 variables set restartfile null

alf/default_scripts/blade_prod.inp

@@ -1,6 +1,13 @@
 ! BLaDE input file for lambda dynamics
 !
 
+! This script expects nsteps (production steps) and itt (chunk of
+! production) to be set by the calling process in arguments.inp
+!
+! Other important variables such as temp (the system temperature) and
+! sysname (the system name in the prep directory) are set by
+! variablesflat.inp
+
 verbose 0
 
 stream arguments.inp

alf/default_scripts/bladelib_flat.inp

@@ -1,6 +1,13 @@
 * CHARMM with BLaDE input file for lambda dynamics
 *
 
+! This script expects esteps (equilibration steps), nsteps (production
+! steps), and seed (random number seed) to be set by the calling process
+!
+! Other important variables such as minimizeflag (whether or not to
+! minimize), temp (the system temperature), and sysname (the system name
+! in the prep directory) are set by variablesflat.inp
+
 set fnex = 5.5
 
 stream "variablesflat.inp"
@@ -20,7 +27,7 @@ endif
 
 faster on
 
-if @minimize .eq. 1 then ! only defined in variables1.inp
+if @minimizeflag .eq. 1 then ! only defined in variables1.inp
 define uninitialized select .not. initialized end
 if ?nsel .eq. 0 then
    ! domdec gpu only dlb off ndir 1 @nnodes 1