plmd

#!/bin/bash
#########################################################################
# plmd, Version 1.5.1
# 2021, ShanghaiTech University, Wang Lin
# Email: wanglin3@shanghaitech.edu.cn
#########################################################################

function checkparam(){
testlogic=false
for line in $2
do
if [ "$1"x == "$line"x ];then
  testlogic=true;break;fi
done
if [ "${testlogic}"x == "true"x ];then
  echo "Note: Parameter "$1" is vaild.";
else
  echo "Error: your parameter "$1" is unvaild, please check it." ;echo ""
  exit
fi
}

function set_msj(){
cat<<MSJ1 > plmd.msj
task {
  task = "desmond:auto"
  set_family = {
    simulate = {
    checkpt.write_last_step = no
    }
  }
}
MSJ1
}

function gen_cfg(){
interval=`echo "scale=1;${mdtime} / ${frames}"|bc`
cat<<MDCFG > plmd.cfg
annealing = false
backend = {
}
bigger_rclone = false
checkpt = {
    first = 0.0
    interval = 240.06
    name = "\$JOBNAME.cpt"
    write_last_step = true
}
cpu = 1
cutoff_radius = 9.0
elapsed_time = 0.0
energy_group = false
eneseq = {
    first = 0.0
    interval = 1.2
    name = "\$JOBNAME\$[_replica\$REPLICA\$].ene"
}
ensemble = {
barostat = {
    tau = 2.0
}
class = NPT
method = MTK
thermostat = {
    tau = 1.0
}
}
glue = solute
maeff_output = {
    first = 0.0
    interval = 120.0
    name = "\$JOBNAME\$[_replica\$REPLICA\$]-out.cms"
    periodicfix = true
    trjdir = "\$JOBNAME\$[_replica\$REPLICA\$]_trj"
}
meta = false
meta_file = ?
pressure = [1.01325 isotropic ]
randomize_velocity = {
    first = 0.0
    interval = inf
    seed = ${seed}
    temperature = "@*.temperature"
}
$constraintDescription
simbox = {
    first = 0.0
    interval = 1.2
    name = "\$JOBNAME\$[_replica\$REPLICA\$]_simbox.dat"
}
surface_tension = 0.0
taper = false
temperature = [
    [${temperature}.0 0 ]
]
time = ${mdtime}
timestep = [0.002 0.002 0.006 ]
trajectory = {
    center = []
    first = 0.0
    format = dtr
    frames_per_file = 250
    interval = ${interval}
    name = "\$JOBNAME\$[_replica\$REPLICA\$]_trj"
    periodicfix = true
    write_velocity = false
}
MDCFG
}

function add_counterion(){
    neutralizeion=$1
cat<<EOF
  add_counterion = {
       ion = ${neutralizeion}
       number = neutralize_system
  }
EOF
};export -f add_counterion

function box(){
    boxshape=$1
    boxsize=$2
cat<<EOF
  box = {
     shape = ${boxshape}
     size = ${boxsize}
     size_type = buffer
  }
EOF
};export -f box

function add_salt(){
    saltconcentration=$1
    positive_ion=$2
    negative_ion=$3
cat<<EOF
  salt = {
     concentration = ${saltconcentration}
     negative_ion = ${negative_ion}
     positive_ion = ${positive_ion}
  }
EOF
};export -f add_salt

function set_sys_build_msj(){
    # genenate parts to system build for msj 
    SYS_array=(${SYS_pipeline//:/ })
    export solvent=${SYS_array[0]}
    export neutralizeion=${SYS_array[1]}
    export salt_buffer_info=${SYS_array[2]}
    [ ${first_input} == "no" ] || checkparam ${neutralizeion} "${all_ions_list[*]}"
    [ ${first_input} == "no" ] || checkparam ${solvent} "${solventlist[*]}"
    echo "build_geometry {" >> plmd.msj
    add_counterion ${neutralizeion} >> plmd.msj
    box ${boxshape} ${boxsize} >> plmd.msj
    echo "  override_forcefield = ${override_forcefield}" >> plmd.msj
    echo "  rezero_system = ${rezerosys}" >> plmd.msj
    echo "  minimize_volume = ${minisys}" >> plmd.msj
    salt_buffer_list=(${salt_buffer_info//+/ })
    for salt_buffer in ${salt_buffer_list[@]} ;do
        salt_buffer_array=(${salt_buffer//-/ })
        export saltconcentration=${salt_buffer_array[0]}
        export positive_ion=${salt_buffer_array[1]}
        export negative_ion=${salt_buffer_array[2]}
        [ ${first_input} == "no" ] || checkparam ${positive_ion} "${positiveionlist[*]}"
        [ ${first_input} == "no" ] || checkparam ${negative_ion} "${nagativeionlist[*]}"
        add_salt ${saltconcentration} ${positive_ion} ${negative_ion} >> plmd.msj
    done
    echo "  solvent = ${solvent}" >> plmd.msj
    echo "}" >> plmd.msj
}

function assign_opls(){
cat<<MSJ2 >> plmd.msj
assign_forcefield {
  forcefield = ${ffid}
  hydrogen_mass_repartition = ${hydrogen_mass_repartition}
  ${water}
  dir      = "."
  compress = ""
}
MSJ2
}

function assign_viparr(){
    input_to_viparr=$1
    viparr_output=$2
    VMD_ASL=`echo ${LigandASL} | sed 's/chain.name/chain/g' | sed 's/res.ptype/resname/g' | sed 's/,/ /g' | sed 's/res.num/resid/g' | sed 's/AND/and/g' | sed 's/NOT/not/g' | sed 's/OR/or/g'`
    water_type=`echo ${solvent} | tr [A-Z] [a-z]`
    if [ ${water_type} == "spce" ] || [ ${water_type} == "tip3p" ] || [ ${water_type} == "tip4pew" ];then
        monovalent_ff="ions.amber1jc.${water_type}"
    else
        monovalent_ff="ions.amber1lm_iod.all"
    fi
    if [ ${ffid} == "Amber" ];then
        viparr ${input_to_viparr} ${viparr_output} -f aa.amber.ff19SB -f aa.amber.phosaa10 -f aa.amber.ffncaa -f aa.amber.ffptm -f ${monovalent_ff} -f lipid.amber.lipid17 -f ions.amber2ff99.tip3p -f na.amber.bsc1 -f na.amber.tan2018 -f water.${water_type} --ligand-selection "${VMD_ASL}"
    elif [ ${ffid} == "DES-Amber" ];then
        viparr ${input_to_viparr} ${viparr_output} -f aa.DES-Amber -f aa.amber.phosaa10 -f aa.amber.ffncaa -f aa.amber.ffptm -f ${monovalent_ff} -f lipid.amber.lipid17 -f ions.amber2ff99.tip3p -f na.amber.bsc1 -f na.amber.tan2018 -f water.${water_type} --ligand-selection "${VMD_ASL}"
    elif [ ${ffid} == "Charmm" ];then
        [ ${water_type} == "tip3p" ] && water_type="tip3p_charmm"
        viparr ${input_to_viparr} ${viparr_output} -f aa.charmm.c36m -f carb.charmm.c36 -f ethers.charmm.c35 -f ions.charmm36 -f lipid.charmm.c36 -f water.${water_type} -f na.charmm.c36 --ligand-selection "${VMD_ASL}"
    fi
}

function set_md(){
# Generate MSJ file
cat<<MD1 >> plmd.msj
simulate {
    title       = "Brownian Dynamics NVT, T = 10 K, small timesteps, and restraints on solute_heavy_atom, 100 ps"
    annealing   = off
    time        = 100
    timestep    = [0.001 0.001 0.003 ]
    temperature = 10.0
    ensemble = {
        class = "NVT"
        method = "Brownie"
        brownie = {
            delta_max = 0.1
        }
    }
    restrain = {
        atom = "solute_heavy_atom"
        force_constant = 50.0
    }
}
simulate {
    title       = "Brownian Dynamics NVT, T = 10 K, small timesteps, and restraints on user defined sets, ${mintime} ps"
    annealing   = off
    time        = ${mintime}
    timestep    = [0.001 0.001 0.003 ]
    temperature = 10.0
    ensemble = {
        class = "NVT"
        method = "Brownie"
        brownie = {
            delta_max = 0.1
        }
        thermostat = {
            tau = 1.0
        }
    }
    ${constraintDescription}
}
simulate {
    title       = "Langevin small steps NVT, T = 10 K, and restraints on solute heavy atoms, 12 ps"
    annealing   = off
    time        = 12
    timestep    = [0.001 0.001 0.003]
    temperature = 10.0
    restrain    = { atom = solute_heavy_atom force_constant = 50.0 }
    ensemble    = {
        class  = NVT
        method = Langevin
        thermostat.tau = 0.1
    }
    randomize_velocity.interval = 1.0
    eneseq.interval             = 0.3
    trajectory.center           = []
}
simulate {
    title       = "Langevin NPT, T = 10 K, and restraints on solute heavy atoms, 12 ps"
    annealing   = off
    time        = 12
    temperature = 10.0
    restrain    = retain
    ensemble    = {
        class  = NPT
        method = Langevin
        thermostat.tau = 0.1
        barostat  .tau = 50.0
    }
    randomize_velocity.interval = 1.0
    eneseq.interval             = 0.3
    trajectory.center           = []
}
simulate {
    title       = "Langevin NPT, T = 10 K, and restraints on solute heavy atoms, 12 ps"
    time        = 12
    restrain    = retain
    ensemble    = {
        class  = NPT
        method = Langevin
        thermostat.tau = 0.1
        barostat  .tau = 50.0
    }
    randomize_velocity.interval = 1.0
    eneseq.interval             = 0.3
    trajectory.center           = []
}
simulate {
    title       = "Langevin NPT and no restraints, 24ps"
    effect_if   = [["@*.*.annealing"] 'annealing = off temperature = "@*.*.temperature[0][0]"']
    time        = 24
    ensemble    = {
        class  = NPT
        method = Langevin
        thermostat.tau = 0.1
        barostat  .tau = 2.0
    }
    eneseq.interval   = 0.3
    trajectory.center = solute
}
simulate {
    title    = "Final MD and analysis, ${mdtime} ps"
    cfg_file = "plmd.cfg"
    jobname  = "\$MASTERJOBNAME"
    dir      = "."
    compress = ""
}
MD1
if [ "$Analysis"x == "true"x ];then
cat<<MSJAA >> plmd.msj
pl_analysis {
    ligand_asl = "$LigandASL"
    protein_asl = "$ProteinASL"
}
MSJAA
fi
}

#Default
    # MD system build
    export SYS_BUILD="yes"
    export System_Mode="INC"
    export negaion='Cl'
    export solvent="SPC"
    export rezerosys=true
    export minisys=true
    export boxshape=cubic
    export boxsize=15.0
    export position=Na
    export ffid=OPLS_2005
    export NeorPo=false
    export saltconcentration=0.15
    export hydrogen_mass_repartition=false
    
    # MD params
    export mintime=100
    export mdtime=100000.0
    export num=1
    export temperature=310
    export frames=1000
    export Analysis=true
    export trunoff_md=false
    export ProteinASL="protein"
    export LigandASL="res.ptype UNK"
    export constraint=false
    export constraintASL="chain.name A AND backbone"
    export constraintForce=10

    # File and Host
    export input=*.mae
    export CPU_HOST=CPU
    export GPU_HOST=GPU

    # Desmond
    export Desmond=${Desmond}
    export viparr=${viparr}
    export VIPARR_FFPATH=${VIPARR_FFPATH}

function forcefield_info(){
cat<<HELP
The current force fields support in plmd:
    S-OPLS: 
        The force fields in Schrödinger packages, recommended to ligand-protein complex.
    OPLS_2005: 
        The default force field of Desmond package.
    Amber: 
        Recommended to protein, DNA, RNA, lipid and other systems.
        Amber-ff19SB for protein, Amber-ffncaa for non-canonical aa, Amber-ffptm for 
        post-translational modifications, amber1jc ion parameters adapt with spce, 
        tip3p or tip4pew, Amber-bsc1 for DNA, Amber-tan2018 for RNA.
    Charmm:
        Recommended to protein, DNA, RNA, lipid, carbohydrate and other systems.
        Charmm36m for protein, Charmm36 for carbohydrate, ions, lipid and nucleic acid.
    DES-Amber:
        Recommended to protein-protein complex.
        DES-Amber for protein-protein complex.
HELP
}

function help(){
cat<<HELP

Usage: plmd [OPTION] <parameter>

An automatic Desmond MD pipline for protein-ligand complex MD simulation.

Example: 
1) plmd -i "*.mae" -S INC -P "chain.name A" -L "res.ptype UNK" -H CPU -G GPU
2) plmd -i "*.mae" -S OUC -P "chain.name A" -L "chain.name B" -t 200 -H CPU -G GPU
3) plmd -i "*.mae" -S "TIP4P:Cl:0.15-Na-Cl+0.02-Fe2-Cl+0.02-Mg2-Cl" -L "res.num 999" -G GPU
4) plmd -i "*.cms" -P "chain.name A" -L "res.ptype ADP" -H CPU -G GPU

Input parameter:
  -i	Use a file name (Multiple files are wrapped in "", and split by ' ') *.mae or *.cms ;
            or regular expression to represent your input file, default is *.mae.

System Builder parameter:
  -S    System Build Mode: <INC>
            INC: System in cell, salt buffer is 0.15M KCl, water is SPC. Add K to neutralize system.
            OUC: System out of cell, salt buffer is 0.15M NaCl, water is SPC. Add Na to neutralize system.
            Custom Instruct: Such as: "TIP4P:Cl:0.15-Na-Cl+0.02-Fe2-Cl+0.02-Mg2-Cl"
                Interactive addition of salt. Add Cl to neutralize system.
                    for positive_ion: Na, Li, K, Rb, Cs, Fe2, Fe3, Mg2, Ca2, Zn2 are predefined.
                    for nagative_ion: F, Cl, Br, I are predefined.
                    for water: SPC, TIP3P, TIP4P, TIP5P, DMSO, METHANOL are predefined.

  -b	Define a boxshape for your systems. <${boxshape}>
            box types: dodecahedron_hexagon, cubic, orthorhombic, triclinic
  -s	Define a boxsize for your systems.  <${boxsize}>
		for dodecahedron_hexagon and cubic, defulat is 15.0;
		for orthorhombic or triclinic box, defulat is [15.0 15.0 15.0];
		If you want use Orthorhombic or Triclinic box, your parameter should be like "15.0 15.0 15.0"
  -R    Redistribute the mass of heavy atoms to bonded hydrogen atoms to slow-down high frequency motions.
  -F	Define a force field to build your systems. <OPLS_2005>
		    OPLS_2005, S-OPLS are recommended to receptor-ligand systems.
            Amber, Charmm, DES-Amber are recommended to other systems. Use -O to show more details.

Simulation control parameter:
  -m	Enter the maximum simulation time for the Brownian motion simulation, in ps. <${mintime}>
  -t    Enter the Molecular dynamics simulation time for the product simulation, in ns. <100>
  -T    Specify the temperature to be used, in kelvin. <${temperature}>
  -N    Number of Repeat simulation with different random numbers. <${num}>
  -P    Define a ASL to receptor, such as "protein".
  -L    Define a ASL to ligand, such as "res.ptype UNK".
  -q    Turn off receptor-ligand analysis.
  -u    Turn off md simulation, only system build.
  -C    Set constraint to an ASL, such as "chain.name A AND backbone"
  -f    Set constraint force, default is 10.
  -o    Specify the approximate number of frames in the trajectory.  <1000>
        This value is coupled with the recording interval for the trajectory and the simulation time: the number of frames times the trajectory recording interval is the total simulation time.
        If you adjust the number of frames, the recording interval will be modified.

Job control:
  -G	HOST of GPU queue, default is ${GPU_HOST}.
  -H    HOST of CPU queue, default is ${CPU_HOST}.
  -D	Your Desmond path. <${Desmond}>
  -V    Your viparr path. <${viparr}>
  -v    Your viparr force fields path. <${VIPARR_FFPATH}>

Thank you for your using, If you found any question, Please contact wanglin3@shanghaitech.edu.cn.
HELP
}

while getopts ":hi:S:nb:s:RF:m:t:T:N:P:L:C:f:o:H:G:D:quV:v:O" opt
do
  case $opt in
    h)
        help
        exit;;
    O)
        forcefield_info
        exit;;
    i)
        export input=$OPTARG;;
    S)
        export System_Mode=$OPTARG;;
    s)
        export boxsize=$OPTARG
        export boxxyz=${boxsize}
        ;;
    b)
        export boxshape=$OPTARG;;
    m)
        export mintime=$OPTARG;;
    R)
        export hydrogen_mass_repartition=true;;
    H)
        export CPU_HOST=$OPTARG;;
    G)
        export GPU_HOST=$OPTARG;;
    D)
        export Desmond=$OPTARG;;
    c)
        export saltconcentration=$OPTARG;;
    f)
        export constraint=true
        export constraintForce=$OPTARG
        ;;
    C)
        export constraint=true
        export constraintASL=$OPTARG
        ;;
    t)
        export mdtime=`echo "${OPTARG} * 1000.0"|bc`;;
    N)
        export num=$OPTARG;;
    q)
        export Analysis=false;;
    P)
        export ProteinASL=$OPTARG
        ;;
    L)
        export LigandASL=$OPTARG
        ;;
    T)
        export temperature=$OPTARG;;
    o)
        export frames=$OPTARG;;
    F)
        export ffid=$OPTARG;;
    u)
        export trunoff_md=true;;
    V)
        export viparr=$OPTARG;;
    v)
        export VIPARR_FFPATH=$OPTARG;;
    ?)
      echo ""
      echo "Error: Do not use undefined options. Try -h."
      echo ""
      exit;;
  esac
done

# Check Desmond and hostfile
if [ -d ${Desmond}/ ];then
  if [ "${Desmond}" == "" ];then
    echo "The environment variable \${Desmond} not found. Please check the Desmond path in your ~/.bashrc."
    exit
  fi
  if [ "`grep -c $CPU_HOST ${Desmond}/schrodinger.hosts`" == "0" ];then
    echo $CPU_HOST "Hostfile not found."
    exit
  fi
  if [ "`grep -c $GPU_HOST ${Desmond}/schrodinger.hosts`" == "0" ];then
    echo $GPU_HOST "Hostfile not found."
    exit
  fi
else
  echo "Desmond not found. Please check your Desmond Path."
  exit
fi

# box type and size
if [ "${boxshape}"x == "orthorhombic"x ];then
  export boxsize="[${boxxyz} ]"
elif [ "${boxshape}"x == "triclinic"x ];then
  export boxsize="[${boxxyz} 60.0 60.0 60.0 ]"
else
  export boxsize=${boxsize}
fi

# MD basic control
seedlist=()
if [ ${num} -gt 3 ];then
    read -p "Do you really want to run more than three random simulations? Make sure you have adequate resources. Press Enter to continue and enter NO to exit!" makesure
    if [ "$makesure" == "NO" ];then
        echo "This is really a wise choice!"
        exit
    else
        echo "Your system will perform ${num} simulations at different random speeds."
        for i in $(seq 1 ${num});do
        seedlist+=("$RANDOM");done
    fi
elif [ "${num}" == "1" ];then
        seedlist=(2007)
        echo "Your system will only perform one simulation."
else
    echo "Your system will perform ${num} simulations at different random speeds."
    for i in $(seq 1 ${num});do
    seedlist+=("$RANDOM");done
fi

################################
## Check and Parse Params
################################
export fflist=("OPLS_2005" "OPLS3e" "OPLS4" "OPLS3" "OPLS2" "Amber" "Charmm" "DES-Amber" "S-OPLS")
export nagativeionlist=("F" "Cl" "Br" "I")
export positiveionlist=("Na" "K" "Li" "Rb" "Cs" "Mg2" "Ca2" "Zn2" "Fe2" "Fe3")
export all_ions_list=("F" "Cl" "Br" "I" "Na" "K" "Li" "Rb" "Cs" "Mg2" "Ca2" "Zn2" "Fe2" "Fe3")
export solventlist=("SPC" "SPCE" "TIP3P" "TIP4P" "TIP5P" "DMSO" "METHANOL")
export boxlist=("orthorhombic" "triclinic" "dodecahedron_hexagon" "cubic")
export waterlist=("SPC" "TIP3P" "TIP4P" "TIP5P")
export first_input="yes"
[ ${first_input} == "no" ] || checkparam $ffid "${fflist[*]}"
[ ${first_input} == "no" ] || checkparam $boxshape "${boxlist[*]}"
if [ "$constraint"x == "true"x ];then
    export constraintDescription="restrain = {atom = \"asl:${constraintASL}\" force_constant = ${constraintForce}}"
else
    export constraintDescription="restrain = none"
fi
for line in ${waterlist[@]}; do
if [ "$solvent" == "$line" ];then
  export water="water = ${solvent}"
  break
fi; done

################################
## Assign Force Field
################################
if [ ${ffid:0:4} == "OPLS" ] || [ ${ffid} == "S-OPLS" ]; then
    export override_forcefield=${ffid}
    VIPPAR_ASSIGN="no"
else
    export override_forcefield="S-OPLS"
    VIPPAR_ASSIGN="yes"
    export PATH=${viparr}:$PATH
    source ${viparr}/activate
fi

################################
## Report MD Information
################################
cat<<OUT

Jobs will run at ${GPU_HOST};

Your MD information is as follows:

MIN time: ${mintime} ps
MD time: ${mdtime} ps
temperature: ${temperature} K
Repeat: ${num}
Rondom numbers list: ${seedlist[@]}

OUT

################################
## Check input file
################################
if ls $input >/dev/null 2>&1;then
  echo "These mae/cms files will be processed:"
  for i in `ls $input`; do echo $i;done
else
  echo "Error: mae/cms file not found, please check it."
  exit
fi
echo ""

for i in `ls $input`; do
    ################################
    ## Check input file and assign system build
    ################################
    if [ ${i##*.} == "cms" ]; then
        SYS_BUILD="no"
    elif [ ${i##*.} == "mae" ];then
        SYS_BUILD="yes"
    else
        echo "ERROR: ${i##*.} is unvaild."
        exit
    fi
    ################################
    ## System Build
    ################################
    input_file_path=`readlink -f $i`
    if [ ${SYS_BUILD} == "yes" ];then
        mkdir -p ${i%%.mae}-SYS_BUILD
        cd ${i%%.mae}-SYS_BUILD
        if [ ${System_Mode} == "INC" ];then
            export SYS_pipeline="SPC:K:0.15-K-Cl"
        elif [ ${System_Mode} == "OUC" ];then
            export SYS_pipeline="SPC:Na:0.15-Na-Cl"
        else
            export SYS_pipeline=${System_Mode}
        fi
        set_msj
        set_sys_build_msj
        ${Desmond}/utilities/multisim -JOBNAME ${i%%.mae}-SYS_BUILD -HOST "${CPU_HOST}:3" -SUBHOST ${GPU_HOST} -maxjob 0 -m plmd.msj -c plmd.cfg -WAIT -lic DESMOND_ACADEMIC:16 -o ../${i%%.*}.cms ${input_file_path}
        cd ..
        system_file=`readlink -f ${i%%.*}.cms`
    fi
    system_file=${input_file_path}
    ################################
    ## Assign Force Fields
    ################################
    if [ ${VIPPAR_ASSIGN} == "yes" ];then 
        assign_viparr ${system_file} "${i%%.*}_viparr.cms"
        if [ ! -f "${i%%.*}_viparr.cms" ];then
            echo "ERROR: Viparr Failed!"
            exit
        fi
        ${Desmond}/run gen_fsys.py "${i%%.*}_viparr.cms"
        system_file=`readlink -f ${i%%.*}_viparr.cms`
    fi
    ################################
    ## Submit MD Simulations
    ################################
    if [ $trunoff_md == true ]; then
        echo "Note: You turn off the MD simulation stage!"
        continue
    fi
    ## Parse Seeds and mkdir
    for seed in ${seedlist[@]};do
        if [ -d ${i%%.*}-${seed}-md ];then
            rm -r ${i%%.*}-${seed}-md
            echo "Remove all old files in ${i%%.*}-${seed}-md."
            mkdir -p ${i%%.*}-${seed}-md
        else
            mkdir -p ${i%%.*}-${seed}-md
        fi
        cd ${i%%.*}-${seed}-md
        gen_cfg # Parse Seeds, temperature and md time
        set_msj # Make a .msj file
        if [ ${VIPPAR_ASSIGN} == "yes" ];then 
            set_md
        else
            assign_opls
            set_md
        fi
        if [ ${GPU_HOST} == "localhost" ];then
            ${Desmond}/utilities/multisim -JOBNAME ${i%%.*}-${seed}-md -HOST "${CPU_HOST}:3" -SUBHOST ${GPU_HOST} -maxjob 0 -m plmd.msj -c plmd.cfg -WAIT -lic DESMOND_ACADEMIC:16 -o ../${i%%.*}-${seed}-md-out.cms ${system_file}
        else
            ${Desmond}/utilities/multisim -JOBNAME ${i%%.*}-${seed}-md -HOST "${CPU_HOST}:3" -SUBHOST ${GPU_HOST} -maxjob 0 -m plmd.msj -c plmd.cfg -lic DESMOND_ACADEMIC:16 -o ../${i%%.*}-${seed}-md-out.cms ${system_file}
        fi
        extend_to_time=`echo "${mdtime} * 2.0"|bc`
cat<<EXTEND > extend.sh
${Desmond}/desmond -JOBNAME "${i%%.*}-${seed}-md" -HOST ${GPU_HOST} -restore ${i%%.*}-${seed}-md.cpt -in ${i%%.*}-${seed}-md-in.cms -cfg mdsim.last_time=${extend_to_time}
EXTEND
cat<<RESUME > resume.sh
${Desmond}/desmond -JOBNAME "${i%%.*}-${seed}-md" -HOST ${GPU_HOST} -restore ${i%%.*}-${seed}-md.cpt -in ${i%%.*}-${seed}-md-in.cms
RESUME
cat<<CLUSTER > cluster.sh
${Desmond}/run trj_cluster.py ${i%%.*}-${seed}-md-out.cms ${i%%.*}-${seed}-md_trj ${i%%.*}-${seed}_cluster -rmsd-asl "$1" -fit-asl "\$2" -n "\$3" -HOST ${CPU_HOST} -JOBNAME ${i%%.*}-${seed}_cluster
CLUSTER
cat<<ANALY > occupancy.sh
${Desmond}/run trj_occupancy.py ${i%%.*}-${seed}-md-out.cms ${i%%.*}-${seed}-md_trj ${i%%.*}-${seed}_occupancy -fit-asl "\${2}" -map-asl "\${1}" -grid-length 50 -JOBNAME ${i%%.*}-${seed}_occupancy
ANALY
cat<<ANALY > ppi.sh
${Desmond}/run analyze_trajectory_ppi.py ${i%%.*}-${seed}-md-out.cms ${i%%.*}-${seed}_ppi "\${1}" "\${2}"
ANALY
        cd ..
        date
        echo "MD system ${i%%.*} with seed ${seed} had been submitted to $GPU_HOST."
    done
done

cat<<FF

All jobs are running! 

NOTE: You can view the *-md.eaf in the job directorys to analysis the trajectories.
NOTE: When the jobs finished, you can run the bash scripts in the job directorys to ....
    Resume previous MD: bash resume.sh 
    Extend MD: bash extend.sh # The Default time to extend is 2-fold of initial time.
    Cluster trajectory: bash cluster.sh  "<rmsd selection>"  "<fit selection>"  "<number>"
    Analysis Occupancy: bash occupancy.sh  "<selection to analysis>"  "<fit selection>"
    Analysis ppi contact: bash ppi.sh  "<Component A>"  "<Component B>" 

Thank you for your using, If you found any question, Please contact wanglin3@shanghaitech.edu.cn.

FF