Allow tinkerio to read TINKERPATH variable from environment

src/liquid.py

@@ -255,6 +255,7 @@ def read_data(self):
  # Read the 'data.csv' file. The file should contain guidelines.
  with open(os.path.join(self.tgtdir,'data.csv'),'rU') as f: R0 = list(csv.reader(f))
  # All comments are erased.
+ R0=[i for i in R0 if len(i)>0] # remove lines without any content first, otherwise it crashes
  R1 = [[sub('#.*$','',word) for word in line] for line in R0 if len(line[0]) > 0 and line[0][0] != "#"]
  # All empty lines are deleted and words are converted to lowercase.
  R = [[wrd.lower() for wrd in line] for line in R1 if any([len(wrd) for wrd in line]) > 0]
@@ -386,8 +387,13 @@ def npt_simulation(self, temperature, pressure, simnum):
  mol2_send = self.mol2
  else:
  mol2_send = []
+ inputs=self.nptfiles + self.scripts + mol2_send
+ #for input in inputs:
+ # if '.xyz' in input:
+ # output=input.replace('.xyz','.arc')
+ # self.extra_output.append(output)
  queue_up(wq, command = cmdstr+' > npt.out 2>&1 ',
- input_files = self.nptfiles + self.scripts + mol2_send + ['forcebalance.p'],
+ input_files = inputs + ['forcebalance.p'],
  output_files = ['npt_result.p', 'npt.out'] + self.extra_output, tgt=self)
 
  def nvt_simulation(self, temperature):
@@ -408,8 +414,14 @@ def nvt_simulation(self, temperature):
  mol2_send = self.mol2
  else:
  mol2_send = []
+ inputs=self.nvtfiles + self.scripts + mol2_send
+ #for input in inputs:
+ # if '.xyz' in input:
+ # output=input.replace('.xyz','.arc')
+ # self.extra_output.append(output)
+
  queue_up(wq, command = cmdstr+' > nvt.out 2>&1 ',
- input_files = self.nvtfiles + self.scripts + mol2_send + ['forcebalance.p'],
+ input_files = inputs + ['forcebalance.p'],
  output_files = ['nvt_result.p', 'nvt.out'] + self.extra_output, tgt=self)
 
  def polarization_correction(self,mvals):

src/nifty.py

@@ -873,6 +873,8 @@ def load_gz():
 
 # Global variable containing a mapping from target names to Work Queue task IDs
 WQIDS = defaultdict(list)
+taskidtooutputfilepathslist={}
+
 
 def getWorkQueue():
  global WORK_QUEUE
@@ -911,19 +913,27 @@ def queue_up(wq, command, input_files, output_files, tag=None, tgt=None, verbose
  @param[in] input_files (list of files) A list of locations of the input files.
  @param[in] output_files (list of files) A list of locations of the output files.
  """
+ outputfilepaths=[]
  global WQIDS
+ global taskidtooutputfilepathslist
  task = work_queue.Task(command)
  cwd = os.getcwd()
  for f in input_files:
  lf = os.path.join(cwd,f)
- task.specify_input_file(lf,f,cache=False)
+ head,f=os.path.split(f)
+ task.specify_file(lf,f, work_queue.WORK_QUEUE_INPUT, cache=False)
  for f in output_files:
  lf = os.path.join(cwd,f)
- task.specify_output_file(lf,f,cache=False)
+ outputfilepaths.append(lf)
+ task.specify_file(f,f, work_queue.WORK_QUEUE_OUTPUT, cache=False)
+
  if tag is None: tag = command
  task.specify_tag(tag)
  task.print_time = print_time
+ if 'GPUDYNAMICS' in os.environ.keys() and 'gas' not in command:
+ task.specify_gpus(1) 
  taskid = wq.submit(task)
+ taskidtooutputfilepathslist[taskid]=outputfilepaths
  if verbose:
  logger.info("Submitting command '%s' to the Work Queue, %staskid %i\n" % (command, "tag %s, " % tag if tag != command else "", taskid))
  if tgt is not None:
@@ -965,6 +975,7 @@ def queue_up_src_dest(wq, command, input_files, output_files, tag=None, tgt=None
 def wq_wait1(wq, wait_time=10, wait_intvl=1, print_time=60, verbose=False):
  """ This function waits ten seconds to see if a task in the Work Queue has finished. """
  global WQIDS
+ global taskidtooutputfilepathslist
  if verbose: logger.info('---\n')
  if wait_intvl >= wait_time:
  wait_time = wait_intvl
@@ -996,6 +1007,7 @@ def wq_wait1(wq, wait_time=10, wait_intvl=1, print_time=60, verbose=False):
  logger.warning("Task '%s' (task %i) failed on host %s (%i seconds), resubmitted: taskid %i\n" % (task.tag, oldid, oldhost, exectime, taskid))
  WQIDS[tgtname].append(taskid)
  else:
+ taskid=task.id
  if hasattr(task, 'print_time'):
  print_time = task.print_time
  if exectime > print_time: # Assume that we're only interested in printing jobs that last longer than a minute.
@@ -1004,6 +1016,24 @@ def wq_wait1(wq, wait_time=10, wait_intvl=1, print_time=60, verbose=False):
  if task.id in WQIDS[tnm]:
  WQIDS[tnm].remove(task.id)
  del task
+ newoutputfilepaths=[] # sort by largest file size and move largest first (like move arc first so dont submit next job before arc and dyn are returned
+ if taskid in taskidtooutputfilepathslist.keys():
+ outputfilepaths=taskidtooutputfilepathslist[taskid]
+ if outputfilepaths!=None:
+ for outputfilepath in outputfilepaths:
+ head,tail=os.path.split(outputfilepath)
+ if os.path.isdir(head):
+ if os.path.isfile(tail):
+ newoutputfilepaths.append(outputfilepath) 
+ newfiles=[os.path.split(i)[1] for i in newoutputfilepaths]
+ newfilestofilepaths=dict(zip(newfiles,newoutputfilepaths))
+ filesizes=[os.stat(thefile).st_size for thefile in newfiles]
+ newfilestofilesizes=dict(zip(newfiles,filesizes))
+ sorteddic={k: v for k, v in sorted(newfilestofilesizes.items(), key=lambda item: item[1],reverse=True)} 
+ for filename in sorteddic.keys():
+ filepath=newfilestofilepaths[filename]
+ shutil.move(os.path.join(os.getcwd(),filename),filepath)
+
 
  # LPW 2018-09-10 Updated to use stats fields from CCTools 6.2.10
  # Please upgrade CCTools version if errors are encountered during runtime.
@@ -1559,7 +1589,7 @@ def process_err(read):
  return Out
 _exec.returncode = None
 
-def warn_press_key(warning, timeout=10):
+def warn_press_key(warning, timeout=1):
  logger.warning(warning + '\n')
  if sys.stdin.isatty():
  logger.warning("\x1b[1;91mPress Enter or wait %i seconds (I assume no responsibility for what happens after this!)\x1b[0m\n" % timeout)

src/tinkerio.py

@@ -314,14 +314,18 @@ def __init__(self, name="tinker", **kwargs):
  self.warn_vn = False
  super(TINKER,self).__init__(name=name, **kwargs)
 
+ def CheckEnvironmentTinker(self):
+ tinkerpath = os.getenv("TINKERPATH", default="")
+ return tinkerpath
+
  def setopts(self, **kwargs):
 
  """ Called by __init__ ; Set TINKER-specific options. """
 
  ## The directory containing TINKER executables (e.g. dynamic)
  if 'tinkerpath' in kwargs:
  self.tinkerpath = kwargs['tinkerpath']
- if not os.path.exists(os.path.join(self.tinkerpath,"dynamic")):
+ if not os.path.exists(os.path.join(self.tinkerpath,"dynamic")) and 'dynamic' not in os.environ.keys():
  warn_press_key("The 'dynamic' executable indicated by %s doesn't exist! (Check tinkerpath)" \
  % os.path.join(self.tinkerpath,"dynamic"))
  else:
@@ -350,8 +354,24 @@ def calltinker(self, command, stdin=None, print_to_screen=False, print_command=F
  # Sometimes the engine changes dirs and the key goes missing, so we link it.
  if "%s.key" % self.name in csplit and not os.path.exists("%s.key" % self.name):
  LinkFile(self.abskey, "%s.key" % self.name)
- prog = os.path.join(self.tinkerpath, csplit[0])
- csplit[0] = prog
+ tinkpath=self.CheckEnvironmentTinker()
+ if tinkpath!=None: # need to override input masterhost tinkerpath, otherwise if worker node has different )S/ libraries then will crash
+ tinkerpath=tinkpath
+ else:
+ tinkerpath=self.tinkerpath
+
+ if 'analyze' in command:
+ gpuvergood=False
+ if stdin=='E,M':
+ gpuvergood=True
+ if command[-1]=='E' or command[-1]=='M':
+ gpuvergood=True
+ csplit[0]='analyze'
+ command=' '.join(csplit)
+ if '_gpu' not in command:
+ prog = os.path.join(tinkerpath, csplit[0])
+ csplit[0] = prog
+
  o = _exec(' '.join(csplit), stdin=stdin, print_to_screen=print_to_screen, print_command=print_command, rbytes=1024, **kwargs)
  # Determine the TINKER version number.
  for line in o[:10]:
@@ -360,12 +380,33 @@ def calltinker(self, command, stdin=None, print_to_screen=False, print_command=F
  if len(vw.split('.')) <= 2:
  vn = float(vw)
  else:
- ls = vw.split('.')
- last = ls[1:]
- first = ls[0]
- vn = first+'.'+''.join(last)
+ ls=vw.split('.')
+ last=ls[1:]
+ first=ls[0]
+ vn=first+'.'+''.join(last)
+ versionstr=''
+ for e in vn:
+ if e.isdigit() or e=='.':
+ versionstr+=e
+ count=versionstr.count('.')
+ if count==2:
+ newversionstr=''
+ for eidx in range(len(versionstr)):
+ e=versionstr[eidx]
+ if eidx>1:
+ if e!='.':
+ newversionstr+=e
+ else:
+ newversionstr+=e
+ versionstr=newversionstr
+ vn=versionstr
+ vn = float(vn)
+ if '_gpu' in command: # then tinker 9 version
+ vn_need=1
+ else:
+ vn_need = 6.3
+
  vn = float(vn)
- vn_need = 6.3
  try:
  if vn < vn_need:
  if self.warn_vn: 
@@ -643,7 +684,7 @@ def evaluate_(self, xyzin, force=False, dipole=False):
  Result = OrderedDict()
  # If we want the dipoles (or just energies), analyze is the way to go.
  if dipole or (not force):
- oanl = self.calltinker("analyze %s -k %s" % (xyzin, self.name), stdin="G,E,M", print_to_screen=False)
+ oanl = self.calltinker("analyze %s -k %s" % (xyzin, self.name), stdin="E,M", print_to_screen=False)
  # Read potential energy and dipole from file.
  eanl = []
  dip = []
@@ -729,6 +770,9 @@ def energy_dipole(self):
  x = "%s.xyz" % self.name
  self.mol.write(x, ftype="tinker")
  Result = self.evaluate_(x, dipole=True)
+ print('Result',Result,flush=True)
+ print('Result["Energy"]',Result["Energy"],flush=True)
+ print('Result["Dipole"]',Result["Dipole"],flush=True)
  return np.hstack((Result["Energy"].reshape(-1,1), Result["Dipole"]))
 
  def normal_modes(self, shot=0, optimize=True):
@@ -910,13 +954,13 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
  if temperature is not None:
  md_defs["integrator"] = "stochastic"
  else:
- md_defs["integrator"] = "beeman"
+ md_defs["integrator"] = "respa"
  md_opts["thermostat"] = None
  # Periodic boundary conditions.
  if self.pbc:
  md_opts["vdw-correction"] = ''
  if temperature is not None and pressure is not None: 
- md_defs["integrator"] = "beeman"
+ md_defs["integrator"] = "respa"
  md_defs["thermostat"] = "bussi"
  md_defs["barostat"] = "montecarlo"
  if anisotropic:
@@ -933,7 +977,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
 
  eq_opts = deepcopy(md_opts)
  if self.pbc and temperature is not None and pressure is not None: 
- eq_opts["integrator"] = "beeman"
+ eq_opts["integrator"] = "respa"
  eq_opts["thermostat"] = "bussi"
  eq_opts["barostat"] = "berendsen"
 
@@ -942,27 +986,31 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
  self.optimize(0, method="bfgs", crit=1)
  os.system("mv %s.xyz_2 %s.xyz" % (self.name, self.name))
  if verbose: logger.info("Done\n")
+ if 'GPUDYNAMICS' in os.environ.keys() and 'gas' not in self.name:
+ dynamickeyword='dynamic_gpu'
+ suffix=''
+ else:
+ dynamickeyword='dynamic'
+ suffix=''
 
  # Run equilibration.
  if nequil > 0:
  write_key("%s-eq.key" % self.name, eq_opts, "%s.key" % self.name, md_defs)
  if verbose: printcool("Running equilibration dynamics", color=0)
  if self.pbc and pressure is not None:
- self.calltinker("dynamic %s -k %s-eq %i %f %f 4 %f %f" % (self.name, self.name, nequil, timestep, float(nsave*timestep)/1000, 
- temperature, pressure), print_to_screen=verbose)
+ self.calltinker(dynamickeyword+" %s -k %s-eq %i %f %f 4 %f %f%s" % (self.name, self.name, nequil, timestep, float(nsave*timestep)/1000,temperature, pressure,suffix), print_to_screen=verbose)
  else:
- self.calltinker("dynamic %s -k %s-eq %i %f %f 2 %f" % (self.name, self.name, nequil, timestep, float(nsave*timestep)/1000,
- temperature), print_to_screen=verbose)
+ self.calltinker(dynamickeyword+" %s -k %s-eq %i %f %f 2 %f" % (self.name, self.name, nequil, timestep, float(nsave*timestep)/1000,temperature), print_to_screen=verbose)
  os.system("rm -f %s.arc" % (self.name))
 
  # Run production.
  if verbose: printcool("Running production dynamics", color=0)
  write_key("%s-md.key" % self.name, md_opts, "%s.key" % self.name, md_defs)
+
  if self.pbc and pressure is not None:
- odyn = self.calltinker("dynamic %s -k %s-md %i %f %f 4 %f %f" % (self.name, self.name, nsteps, timestep, float(nsave*timestep/1000), 
- temperature, pressure), print_to_screen=verbose)
+ odyn = self.calltinker(dynamickeyword+" %s -k %s-md %i %f %f 4 %f %f%s" % (self.name, self.name, nsteps, timestep, float(nsave*timestep/1000),temperature, pressure,suffix), print_to_screen=verbose)
  else:
- odyn = self.calltinker("dynamic %s -k %s-md %i %f %f 2 %f" % (self.name, self.name, nsteps, timestep, float(nsave*timestep/1000), 
+ odyn = self.calltinker(dynamickeyword+" %s -k %s-md %i %f %f 2 %f" % (self.name, self.name, nsteps, timestep, float(nsave*timestep/1000), 
  temperature), print_to_screen=verbose)
 
  # Gather information.
@@ -977,25 +1025,39 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
  edyn.append(float(s[2]))
  if 'Current Kinetic' in line:
  kdyn.append(float(s[2]))
- if len(s) > 0 and s[0] == 'Temperature' and s[2] == 'Kelvin':
- temps.append(float(s[1]))
+ temps.append(temperature) # GPU version does not print out Temperature (or just new dynamic doesnt)
 
  # Potential and kinetic energies converted to kJ/mol.
  edyn = np.array(edyn) * 4.184
  kdyn = np.array(kdyn) * 4.184
  temps = np.array(temps)
 
+ mass = 0.0
  if verbose: logger.info("Post-processing to get the dipole moments\n")
- oanl = self.calltinker("analyze %s-md.arc" % self.name, stdin="G,E,M", print_to_screen=False)
+ oanl = self.calltinker("analyze %s " % self.name, stdin="G,E,M", print_to_screen=False)
+ for ln, line in enumerate(oanl):
+ if 'Polarization' in line:
+ if 'MUTUAL' in line or 'THOLE' in line:
+ continue
+
+ strip = line.strip()
+ s = line.split()
+ if 'Total System Mass' in line:
+ mass = float(s[-1])
+
+ oanl = self.calltinker("analyze %s-md.arc" % self.name, stdin="E,M", print_to_screen=False)
 
  # Read potential energy and dipole from file.
  eanl = []
  dip = []
- mass = 0.0
  ecomp = OrderedDict()
  havekeys = set()
  first_shot = True
  for ln, line in enumerate(oanl):
+ if 'Polarization' in line:
+ if 'MUTUAL' in line or 'THOLE' in line:
+ continue
+
  strip = line.strip()
  s = line.split()
  if 'Total System Mass' in line:
@@ -1007,10 +1069,13 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None,
  if first_shot:
  for key in eckeys:
  if strip.startswith(key):
+
+
+ value=float(s[-2])*4.184
  if key in ecomp:
- ecomp[key].append(float(s[-2])*4.184)
+ ecomp[key].append(value)
  else:
- ecomp[key] = [float(s[-2])*4.184]
+ ecomp[key] = [value]
  if key in havekeys:
  first_shot = False
  havekeys.add(key)