duck_cli.svl

#svl
//	duck_cli.svl : Command Line for setting up a dynamic undocking calculation
//  04 Feb 2017 (PS) Adapted for command line usage 
//	
// COPYRIGHT (C) 1996-2016 CHEMICAL COMPUTING GROUP INC.  ALL RIGHTS RESERVED.
//
// PERMISSION TO USE, COPY, MODIFY AND DISTRIBUTE THIS SOFTWARE IS HEREBY
// GRANTED PROVIDED THAT: (1) UNMODIFIED OR FUNCTIONALLY EQUIVALENT CODE
// DERIVED FROM THIS SOFTWARE MUST CONTAIN THIS NOTICE; (2) ALL CODE DERIVED
// FROM THIS SOFTWARE MUST ACKNOWLEDGE THE AUTHOR(S) AND INSTITUTION(S); (3)
// THE NAMES OF THE AUTHOR(S) AND INSTITUTION(S) NOT BE USED IN ADVERTISING
// OR PUBLICITY PERTAINING TO THE DISTRIBUTION OF THE SOFTWARE WITHOUT
// SPECIFIC, WRITTEN PRIOR PERMISSION; (4) ALL CODE DERIVED FROM THIS SOFTWARE
// BE EXECUTED WITH THE MOLECULAR OPERATING ENVIRONMENT (MOE) LICENSED FROM
// CHEMICAL COMPUTING GROUP INC.
//
// CHEMICAL COMPUTING GROUP INC. DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS
// SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
// AND IN NO EVENT SHALL CHEMICAL COMPUTING GROUP INC. BE LIABLE FOR ANY
// SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
// RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF
// CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
// CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

#set title	'Dynamic UndoCKing'
#set class	'MOE:simulation'
#set version	'2018.11'

#set main	'runDuck' 



global JPATH;
global FFPATH;
global SCRIPTSPATH;


function _Atoms;
function fwrite_MOE;
function fwrite_AmberMOL2;
function export2FRCMod;
function writeAmberInputs;
function launchAmber;
//SRC MINO
function writeAmberLaunchScripts;
function launchAmber_MINO;
function launchAmber_SLURM;
function launchAmber_Workstation;
function writeAmberLaunchScripts_MINO;
function writeAmberLaunchScripts_SLURM;
function my_fwrite_PDB;
function MM;
function PartialCharge;
function ph4_aDonor, ph4_aAcceptor;



// estimate the memory needed for the calculation

// returns all the atoms from set a that are within 'dist' angstrom from b
local function atoms_within[ a, b, dist ]
    if (length b) == 0 or (length a) == 0 then
	return [];
    endif
    local pkey = prox_open [ dist, aPos b, dist ];
    local [seg,idx,sqr_dist] = prox_find[ pkey, aPos a, 0.0 ];
    prox_close pkey;
    local aidx =  x_pack seg;
    return a[aidx];
endfunction





// stolen from prolib.svl
local function replace_string [s, pos, replacement]
    local s1, s2;
    s1 = keep [s, (first pos) -1];
    s2 = keep [s, -(length s - (last pos))];
    return cat [s1, replacement, s2];
endfunction

// stolen from prolib.svl
local function inside_bracket [s,p]
    local bc = strpos ["]", drop [s, p]];
    if bc > 0 then bc = bc + p; endif
    local bo = strpos ["[", drop [s, p]];
    if bo > 0 then bo = bo + p; endif
    local c;
    if bo == 0 and bc <> 0 then
	c = 1;
    elseif bc == 0 and bo <> 0 then
	c = 1;
    else
	c = select [1,0, bc < bo];
    endif
    return c and not (s(p) == "[" or s(p) == "]");
endfunction

// stolen from prolib.svl
local function search_and_replace_all [s, word, replacement, opt]
    opt = cat [opt, [add_brackets:0]];
    local p, rplmnt;
    if word === "" then return s; endif
    while p = strpos [word, s] loop
	if opt.add_brackets === 1 then
	   if not inside_bracket [s, p] then
		rplmnt = cat ["[", replacement, "]"];
	   else
		rplmnt = replacement;
	   endif
	endif
	s = replace_string [s, [p, p + length word - 1 ], rplmnt];
    endloop
    return s;
endfunction


local function makeMDDirs[]
	local fileList=fbase ftail flist [];
	if indexof['lib',fileList]==0 then
		fmkdir 'lib';
	endif
endfunction

local function writeDoLeap[ligand,dirName,options,suppFFParams,lenRec,solvation]
	local fh;
	if(solvation==1) then
	    fh=fopenw twrite['{}/do.leap',dirName];	//open file handle
	else
	    fh=fopenw twrite['{}/doleaped.leap',dirName];	//open file handle
	endif
	print options.ffParams;
	//please leave both options here so in any case we retain compatibility between old amber versions and newer ones
//	fwrite [fh,'source oldff/leaprc.ff99SBildn\n'];
//	fwrite [fh,'source leaprc.ff99SBildn\n'];
	fwrite [fh,'source oldff/leaprc.ff14SB\n'];
	fwrite [fh,'source leaprc.water.tip3p\n'];
//	fwrite [fh,'loadOff ions94.lib\n'];

	if length options.ffParams[1]>0 then
		local el;
		local mkeyList=[];
		local mkey;
		local paramFileName;
		for el in cat options.ffParams[1] loop
			mkeyList=cat [mkeyList,Message[0,twrite['Adding additional parameters to leap file: {}',el]]];
			paramFileName=(cat suppFFParams) | rotr (cat suppFFParams)==el;
			fwrite [fh,'loadOff {}/{}.lib\n',FFPATH,paramFileName];
			fwrite [fh,'loadAmberPrep {}/{}.prp\n',FFPATH,paramFileName];
		endloop
		sleep 2;
		for mkey in mkeyList loop
			Message[mkey,''];
		endloop
	endif

	if length ligand>0 then
//		local resNames=uniq rName oResidues ligand;
		local res;
		for res in uniq oResidues ligand loop
			local ligName=rName res;
			print cat['Identified Residue ',ligName,' as ligand'];
			fwrite [fh,'loadAmberParams {}.frcmod\n',ligName];
			fwrite [fh,'{}=loadMol2 {}_ok.mol2\n',ligName,ligName];
		endloop
		fwrite [fh,'receptor=loadPdb system.pdb\n'];
		fwrite [fh,'system=combine \{{receptor'];

		for res in uniq oResidues ligand loop
			fwrite [fh,' {}',rName res];
		endloop
		fwrite [fh,'\}\n'];
	else 
		fwrite [fh,'system=loadPdb system.pdb\n'];
	endif
	if(solvation==1) then
	  fwrite [fh,'addIons system Cl- 0\n'];
	  fwrite [fh,'addIons system Na+ 0\n'];
	  fwrite [fh,'solvateBox system TIP3PBOX 18.0 0.65\n'];
	  fwrite [fh,'saveAmberParm system system_solv.prmtop system_solv.prmcrd\n'];
	  fwrite [fh,'savePdb system system_solv.pdb\n'];
	else
	  fwrite [fh,'savePdb system system_leaped.pdb\n'];
	  
	endif
	fwrite [fh,'quit'];

	
	fclose fh;	//close file handle
	
endfunction


local function writeSystemPdb[receptor,fName]
	local pdb_options = [
	    use_chain_id    : 1,
	    force_TER	    : 1
	];
	pdb_options.hnames = 'PDB v3.0';
	pdb_options.amber = 1;

	my_fwrite_PDB [fName,receptor,pdb_options];
endfunction

local function writeLigand[ligand,dirName]
	if length ligand == 0 then
		return [];
	endif
	local count=0;
	local res= uniq oResidues ligand;
	local curRes;
	for curRes in res loop
		print "writing";
		local resName=rName curRes;
		

		local selection=cat oAtoms curRes;
		print resName;
		print length selection;
		export2FRCMod[tok_cat[dirName,resName,'.frcmod'],selection];
		fwrite_AmberMOL2[tok_cat[dirName,resName,'_ok.mol2'], selection, []];
	endloop
	//local resName=first rName oResidues ligand;
	
endfunction

local function executeLeap [path,leapfile]
	local currentDir= cd '.';
	cd path;
	local pid=exe_open [twrite['tleap -f {}',leapfile]];
	local fout=exe_stdout pid;
	local ferr=exe_stdout pid;
	local tleap_msg_key = Message [0, 'Configuring the system using Ambertools'];
	loop
		sleep 1;
	until exe_status pid ==0
	endloop
	exe_close pid;
	Message[tleap_msg_key,[]];
	
	local msgKey=Message [0,'Tleap finished running.'];
	
	
	
	
	cd currentDir;
	
	return msgKey;
	
endfunction

local function get_HB_atomKey[res_ID,atom_Name]
        local res_list=uniq cat aResidue [Atoms[]];
        local mask_res=cat rUID [uniq cat aResidue [Atoms[]]] == res_ID;
        local sel_res=res_list | mask_res;
        local mask_atom=cat aName[cat rAtoms[sel_res]] == atom_Name;
        local sel_atoms=cat rAtoms[sel_res];
        local finalatomkey=sel_atoms | mask_atom;
        return finalatomkey;
endfunction


//local function DEPR_selectAtomsFromReceptorAndLigandDB[selAtomPosition, options]
//	local a=Atoms[];
	
// Select the reference atom from coordinates
//loop over all selected atoms if multiple atoms selected
// 	for curPos in selAtomPositions loop
//	    aSetSelected[a[x_min app add tr abs(curPos-aPos Atoms[])],1];
//  	endloop
//	aSetSelected[a[x_min app add tr abs(selAtomPosition-aPos Atoms[])],1];
// XB - ini
//	changed this expression for a simpler one (not sure what this $$ligand thing is...)
//	local ligand_atoms=_Atoms '$$ligand' | (rName aResidue _Atoms '$$ligand' == 'LIG');
// PS - comment, this is some stuff that can be useful in some cases but not for this application here.

// 	sqrt(add(pow [(oCenterOfMass s2 - oCenterOfMass s1),2 ]))

// 	local ligand_atoms = a | rName aResidue a == 'LIG';
//	keep only the complementary atoms of the ligand
  // 	local asubset;
  // 	if isdonor==1 then
  // 		 asubset = ligand_atoms | ph4_aAcceptor ligand_atoms;
  // 	else
  // 		 asubset = ligand_atoms | ph4_aDonor ligand_atoms;
  // 	endif
  // 	if length asubset < 1 then
  // 		print 'This ligand does not have any atom of the right type.';
  // 	endif
  // 	aSetSelected[asubset[x_min aDist [SelectedAtoms[],asubset]],1];
  // 	if add aSelected a <> 2 then
  // 		print 'Number of selected atoms different than 2. Something went wrong. Must end here.';
  // 		exit[];
  // 	endif
  // 	
  // // PS not used anymore
  // //	local dist = aDist [first SelectedAtoms[], second SelectedAtoms[]];
  // //	if dist > 3.2 then
  // //		print 'Selected atoms are not within H-bond distance. Must end here.';
  // //		exit [];
  // //	endif
  // // XB - end
  // 
  // 	
  // endfunction

function startMD;

local function startDB[options,suppFFParams]
	pot_Load '$MOE/lib/pfrosst.ff';


	local currentDir=cd '.';
	
//SRC changes to get atom key of the HBond residue	
	local selected_AtomKey=get_HB_atomKey[options.HB_resID,options.HB_atomName];
        print [selected_AtomKey];
//SRC end
	local selAtomPosition=aPos selected_AtomKey;
	local tmpName='tmp_receptor.moe';
	SaveAs tmpName;
	
	Close[force:1];
	//copy some important launch and analysis scripts 
//SRC moe2016
	JPATH = getenv 'HOME';
        FFPATH = tok_cat [ JPATH , '/moefiles/svl/tleap'];
        SCRIPTSPATH = tok_cat [ JPATH , '/moefiles/svl/scripts'];
	//JPATH = getenv 'MOE_SVL_LOAD';
    	//SCRIPTSPATH = tok_cat [ JPATH , '/scripts'];
//SRC end
	//local rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/duck_template_gpu.sh'], 'duck_template_gpu.sh'];
	//rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/duck_template_325K_gpu.sh'], 'duck_template_325K_gpu.sh'];
	local rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/getWqbValues.R'], 'getWqbValues.R'];

        if options.queue_IN === 'Minotauro' then
        	rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/submit_duck_smd_gpu_MINO.csh'], 'submit_duck_smd_gpu.csh'];
        endif
        if options.queue_IN === 'SLURM' then
        	rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/submit_duck_smd_gpu_SLURM.csh'], 'submit_duck_smd_gpu.csh'];
        endif
        if options.queue_IN === 'Marc' then
                rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/submit_duck_smd_gpu_MARC.csh'], 'submit_duck_smd_gpu.csh'];
        endif
	if options.input_WS == 1 then
		rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/launch_duck_gpu_Workstation.sh'], 'launch_duck_gpu_Workstation.sh'];
	endif
	
	
	
	print 'Opening ligand MDB file';
	if ftype options.ligand_db<>'file' then 
		print 'Opening ligand database not possible. Please check that your path is correctly pointing to a valid mdb file.';
		return;
	endif
	local mdb=db_Open options.ligand_db;
	local entries=db_Entries mdb;	

	local molfield='mol';
	local [fields, types] = db_Fields mdb;
	local molfieldtype = types | (fields == molfield);
	if molfieldtype<>'molecule' then 
		write ['Warning, field named {} does not exist.',molfield];
		molfield=db_FirstFieldType[mdb,'molecule'];
		write [' Using field named {} instead.',molfield];
	endif

	local entry;
	//loop over all entries to get the ligand molecule
	local molCount=1;

//removed for avoiding writing q files in main folder
////SRC MINO queue in
//	//writeAmberLaunchScripts options;
//	print options.queue_IN;
//	if options.queue_IN === 'Minotauro' then
//		writeAmberLaunchScripts_MINO options;
//	endif
// if options.queue_IN === 'SLURM' then
//		writeAmberLaunchScripts_SLURM options;
//	endif
//	if options.queue_IN === 'Marc' then 
//		writeAmberLaunchScripts options;
//	endif
////SRC MINO queue end
	local molCount_mopacError=0; //SRC define for MOPAC
	local fh_MOPAC_created=0; //SRC define for MOPAC
	for entry in entries loop
		Open tmpName;
		local msgKey=Message[0,tok_cat['Preparing simulations for molecule ',totok molCount]];
		local protein=cat oAtoms Chains[];

		local ligand=mol_Create first db_ReadFields [mdb,entry,molfield];

//		local ligandName=first first ligand;

// get rid of minimization
//		if options.tethered_min_DB==1 then
//			//set thether only on protein atoms
//			print 'setting tethers to ligand and receptor';
//			local heavyProteinAtoms=protein | ((aAtomicNumber protein)<>1);
//			aSetTether[heavyProteinAtoms, aPos heavyProteinAtoms,10,0,0.25];
//			local heavyLigandAtoms=cat oAtoms ligand | ((aAtomicNumber cat oAtoms ligand)<>1);
//			aSetTether[heavyLigandAtoms, aPos heavyLigandAtoms,1,0,0.25];
//			print 'starting minimisation';
//			MM[pot_charge:1];
//		endif
		
		//rename residues systematically to LIG
		local ligResidues=uniq cat aResidue oAtoms ligand;
		local el;

	        for el in ligResidues loop
                        rSetName[el,'LIG'];
                        local ligAtoms=cat oAtoms el;


//SRC MOPAC ERROR
                	local flag_mopacError = 0;
			local [result,return_code] = task_call['PartialCharge',[ligAtoms, 'AM1-BCC'],[errmsg:'ignore']];
			if return_code == 'error' then
				flag_mopacError=1;
			else
				local q=first result;
				aSetCharge [ligAtoms,q];
			endif	
                endloop

                if flag_mopacError == 1 then
			print 'Error calculating charges, skipping ligand.';
                        if molCount_mopacError == 0 then
				fh_MOPAC_created=1;
                                local fh_MOPAC;
                                fh_MOPAC = fopenw twrite['DUCK_charge_errors.txt','./'];   //open file handle
                                molCount_mopacError = 1;
                        endif;
                        fwrite [fh_MOPAC,'LIG_target_{}\n',molCount];
                        molCount = molCount+1;
                        Close[force:1];
                        Message[msgKey,[]];
                        continue;
                endif;
//SRC MOPAC ERROR end	
		
		local fileList=fbase ftail flist '.';
		//create folder name for writing simulation parameters
		local ligandDirName=tok_cat['LIG_target_',totok molCount];
		
		if indexof[ligandDirName,fileList]==0 then
			fmkdir ligandDirName;
		endif
		
		cd ligandDirName;
		local receptorSelection=[];
		local a=Atoms[];
		local curPos;

		for curPos in (tr selAtomPosition) loop
		    print curPos-aPos Atoms[];
		    receptorSelection=cat [receptorSelection,a[x_min app add tr abs(curPos-aPos Atoms[])]]; //instead of selecting add it to a set
		endloop
		
		local ligandSelection=[];
		local receptorCOM=oCenterOfMass receptorSelection;	//calculate center of mass of receptor selection
		
		local ligand_atoms = a | rName aResidue a == 'LIG';
		//ligand_atoms=ligand_atoms | aAtomicNumber ligand_atoms <> 1;	//get only heavy atoms
		local m1 = aAtomicNumber ligand_atoms <> 1;
        local m2 = aAtomicNumber ligand_atoms <> 6;
        ligand_atoms=ligand_atoms | m1*m2;      // get only non-hydrogen and non-carbon atom

		local distances=cat sqrt add pow[tr(tr aPos ligand_atoms - [receptorCOM]),2];
		// by default now use closest atom on ligand
		ligandSelection=ligand_atoms[x_min distances];


		// uncomment this if you want to reintroduce the environment options from the moe window

		//if options.environment=='Closest Atom on Ligand' then
		//    ligandSelection=ligand_atoms[x_min distances];
		//else 
		//    ligandSelection=ligand_atoms | distances<=options.smd_cutoff;
		    
		//endif

		
		//selectAtomsFromReceptorAndLigandDB[selAtomPosition,isdonor];	//TODO check this
		startMD[cat oAtoms ligand,options,suppFFParams,ligandSelection,receptorSelection]; // TODO CONTINUE
		cd '../';
		
		
		molCount=molCount+1;
		Close[force:1];
		Message[msgKey,[]];	
		
	endloop
	if fh_MOPAC_created == 1 then
		fclose[fh_MOPAC];
	endif;
	db_Close mdb;
endfunction

local function getSelectedAtomIdsFromPDB[path,atomSelection]
	local maxANumber = max aNumber Atoms[];
	local newChains = ReadAuto [ path, [center:0,select:0], 'pdb' ];
	local newAtoms=cat oAtoms newChains;
	local curAtom;
	local aIndices; 
	local distances;
	local newAtomSelection;
	for curAtom in atomSelection loop
	    distances=aDist [newAtoms,curAtom];
	    newAtomSelection=cat[newAtomSelection,newAtoms[x_min cat distances]];
	
	endloop;
	newAtomSelection=uniq newAtomSelection; 
// 	local selAtom1=first SelectedAtoms[];
// 	local selAtom2=second SelectedAtoms[];
// 
// 	local diff1=aDist [selAtom1,newAtoms];
// 	local diff2=aDist [selAtom2,newAtoms];
// 	
// 	local firstIndex=aNumber newAtoms[x_pack(diff1<0.01)];
// 	local secondIndex=aNumber newAtoms[x_pack(diff2<0.01)];
// 	print firstIndex;
// 	print secondIndex;
// 	print aIndices;
	local ligand_selection=newAtomSelection | rName aResidue newAtomSelection == 'LIG';
	local receptor_selection=newAtomSelection | rName aResidue newAtomSelection <> 'LIG';

	local ligand_com=oCenterOfMass ligand_selection;
	local receptor_com = oCenterOfMass receptor_selection;
	
	local comdist=sqrt add pow [ligand_com-receptor_com,2];
	local ligand_aNumbers=aNumber ligand_selection;
	local receptor_aNumbers=aNumber receptor_selection;
 	oDestroy newChains;
	return([ligand_aNumbers-maxANumber,receptor_aNumbers-maxANumber,comdist]);
endfunction



local function startMD[ligand,options,suppFFParams,ligandSelection,receptorSelection]
	print 'Using Prm@frosst forcefield';
	pot_Load '$MOE/lib/pfrosst.ff';
	local currentDir=cd '.';
	print cat ['Creating folders in ',currentDir];
	makeMDDirs[];
	local allAtoms=Atoms[];
	local receptor=[];

	local ligandAtomNames=aName ligand;
	aSetName [ligand, aMMType ligand];	//overwrite ligand names with MMtype names for mol2 export
	receptor=diff[allAtoms,ligand];	//from all atoms filter out the ligand atoms to retain only resting atoms
	//treat ligand and receptor
	print cat['Selected ',length ligand, ' ligand atoms'];

	print cat['Selected ',length receptor, ' receptor atoms'];	
	writeSystemPdb [receptor,'lib/system.pdb'];
	writeLigand[ligand,'lib/'];
	if length ligand>0 then
		aSetName [ligand, ligandAtomNames];	//rewrite original atom names back to atoms
	endif
	writeDoLeap[ligand,'lib',options,suppFFParams,length receptor,1];
	writeDoLeap[ligand,'lib',options,suppFFParams,length receptor,0];
	local msgKey=executeLeap['lib','do.leap'];
	
	
	//check if successful
	if fsize 'lib/system_solv.prmtop'== 0 then
		Message[msgKey,[]];
		msgKey=Message[0, 'Error, tleap did not manage to build the system'];
		print 'Error, tleap did not manage to build the system';
		sleep 5;

		return [];
	endif
	Message[msgKey,[]];

	//rerun tleap to get proper prmtop atom numberings

	msgKey=executeLeap['lib','doleaped.leap'];
	local selectedAtomIds=getSelectedAtomIdsFromPDB['lib/system_leaped.pdb',cat [ligandSelection,receptorSelection]];
	Message[msgKey,[]];
	local lenRec = length uniq aResidue oAtoms receptor;
	local lenOtherAtoms = app add [ rType uniq aResidue oAtoms receptor == 'none'];
	lenRec=lenRec-lenOtherAtoms;
// 	local selAtomsPositions=tr aPos SelectedAtoms[];
// 	local distSelected= aDist[first SelectedAtoms[], second SelectedAtoms[]];
	local ret=writeAmberInputs['',options,ligand,lenRec,selectedAtomIds];

//SRC COMMENT AS FOLLOWING LOOP DECIDES WHICH LAUNCHAMBER TO CALL
//	if ret==0 then
//		launchAmber[options];
//	else
//		sleep 3;
//		Message[ret,''];
//	endif
/////

//SRC MINO or Marc queues
	if options.queue_IN === 'Minotauro' then
		if ret==0 then
			launchAmber_MINO options;
		else
			sleep 3;
			Message[ret,''];
		endif
	endif
	if options.queue_IN === 'SLURM' then
		if ret==0 then
			launchAmber_SLURM options;
		else
			sleep 3;
			Message[ret,''];
		endif
	endif
	if options.queue_IN === 'Marc' then 
		if ret==0 then
			launchAmber options;
		else
			sleep 3;
			Message[ret,''];
		endif
	endif
	//SRC if local workstation also run
	//print options.input_WS;
	if options.input_WS == 1 then
		local rcode = _fcopy [tok_cat [ SCRIPTSPATH,'/launch_duck_gpu_Workstation.sh'], '../launch_duck_gpu_Workstation.sh'];
		launchAmber_Workstation options;
	endif
//SRC MINO or Marc queues end
endfunction


local function loadSuppFFParams[]
	local fh=fopenr[twrite['{}/register.txt',FFPATH]];
	local tmpText;
	tmpText=fread[fh,'{G w=`\t` r=`\n`}{t:}{t:}'];
	local titles=[];
	if tmpText[2]>0 then
		titles=cat [tmpText[1][1]];
	endif
	while tmpText[2]>0 loop
		tmpText=fread[fh,'{G w=`\t` r=`\n`}{t:}{t:}'];		
		if tmpText[2]>0 then
			titles=cat [titles,tmpText[1][1]];
		endif
	endloop
	fclose fh;
	return titles;
endfunction


//run duck with a tagged vector giving 
//moefile: input receptor file with one selected atom
//mdbfile: input mdb with all ligands
global function runDuckCLI opt
//SRC moe 2016
	JPATH = getenv 'HOME';
        FFPATH = tok_cat [ JPATH , '/moefiles/svl/tleap'];
        SCRIPTSPATH = tok_cat [ JPATH , '/moefiles/svl/scripts'];
    	//JPATH = getenv 'MOE_SVL_LOAD';
    	//FFPATH = tok_cat [ JPATH , '/tleap'];
//SRC end
	Open opt.moefile;

	local suppFFParams=loadSuppFFParams[];
	local md_options=['ffParams':[ [], 0 ], 'queue_IN':opt.queue, 'input_WS':0 ,'ligand_db':opt.mdbfile,'smd_fc':50,'smd_length':0.5,'smd_displ':2.5,'smd_cutoff':5.0,'equilibration_length':1,'md_length':0.5,'wqb_thresh':opt.wqbThreshold,'wqb_maxruns':opt.wqbMaxRuns,'HB_resID':opt.HB_resID,'HB_atomName':opt.HB_atomName];
	startDB[md_options,suppFFParams];		
    
endfunction

//HOW TO EXECUTE:
//moebatch -load duck_cli.svl -exec "runDuckCLI[moefile:'rec.moe',mdbfile:'small.mdb',queue:'Marc', wqbThreshold:6.0,wqbMaxRuns:10,HB_resID:66,HB_atomName:'N']";

#eof