# YASARA MACRO
# TOPIC:       3. Molecular Dynamics
# TITLE:       Running a molecular dynamics simulation of a membrane protein
# REQUIRES:    Dynamics
# AUTHOR:      Elmar Krieger
# LICENSE:     GPL
# DESCRIPTION: This macro sets up and runs a membrane simulation. It scans the protein for secondary structure elements with hydrophobic surface residues, orients it accordingly and embeds it in a membrane of adjustable lipid composition. Finally a 250 ps restrained equilibration simulation is run, which ensures that the membrane can adapt to the newly embedded protein. Then the real simulation starts.

# Parameter section - adjust as needed, but NOTE that some changes only take
# effect if you start an entirely new simulation, not if you continue an existing one. 
# ====================================================================================

# The structure to simulate must be present with a .pdb or .sce extension.
# If a .sce (=YASARA scene) file is present, the membrane and cell must have been added.
# You can either set the target structure by clicking on Options > Macro > Set target,
# or by uncommenting the line below and specifying it directly.
#MacroTarget = 'c:\MyProject\1crn'

# Extension of the cell on each side of the protein in the membrane plane (=XZ plane)
# '15' means that the membrane will be 30 A larger than the protein
memextension=15

# Extension of the cell on each side of the protein along the third (water) axis (=Y-axis)
# '10' means that the cell will be 20 A higher than the protein
waterextension=10

# Membrane composition: The percentages of phosphatidyl-ethanolamine (PEA),
# phosphatidyl-choline (PCH) and phosphatidyl-serine (PSE), must sum up to 100.
# Note that PCH has a large headgroup which cannot form hydrogen bonds, and thus
# reduces membrane stability. PEA is the most stable membrane lipid.
# Percentage of phosphatidyl-ethanolamine
PEApercent=100
# Percentage of phosphatidyl-choline
PCHpercent=0
# Percentage of phosphatidyl-serine
PSEpercent=0
# Or uncomment below to use your own membrane, see membrane simulation docs for details
#usermemname='YourChoice'
#usermemsize=77.21,73.24

# pH at which the simulation should be run, by default physiological pH
ph=7.4

# NaCl concentration in mass percent (0.9% is a physiological solution)
nacl=0.9

# Forcefield to use (this is a YASARA command, so no '=' used)
ForceField AMBER03

# Simulation temperature
temperature='298K'

# Name of solvent molecules used to measure solvent density.
# If you use mixed solvents, check the PressureCtrl documentation. 
solvent='HOH'

# Solvent density in g/ml (0.997 is water at 298K)
# If you run at a significantly different temperature, this density needs to
# be adapted. Look in the PressureCtrl documentation.
# For solvents other than water, you have to create your own solvent box
# as described in the FillCellObj documentation and save it as .._solvent.sce.
density = 0.997

# Pressure at which the simulation should be run [bar].
# This is combined with the solvent density as described in the PressureCtrl docs.
pressure=1 

# Cutoff
cutoff=7.86

# Use PME for longrange electrostatics
longrange='Coulomb'

# Equilibration period in picoseconds:
# During this initial equilibration phase, the membrane is artificially stabilized
# so that it can repack and cover the solute, while solvent molecules are kept outside.
equiperiod=250

# Duration of the complete simulation, must be longer than equiperiod above.
# Alternatively use e.g. duration=5000 to simulate for 5000 picoseconds
duration='forever'

# Delay for animations, 1=maximum speed
delay=100

# Constrain bond lengths to hydrogens and water bond angles to allow a larger timestep
# If set to 'yes', the MD will run faster, but will be a bit less accurate
constrain='no'

# The format used to save the trajectories, sim or xtc
format='sim'

# Normally no change required below this point
# ============================================

RequireVersion 9.9.25

# Keep the solute from diffusing around and crossing periodic boundaries
CorrectDrift On

# Treat all simulation warnings as errors that stop the macro
WarnIsError On

# Membrane simulations are always periodic
Boundary periodic

# The membrane values below should not be changed, since they only describe
# membrane characteristics that originate from the simulation: during the
# simulation, the membrane reaches an equilibrium, which depends only on
# the lipid composition, the embedded protein and the force field parameters,
# but not on the initial setup. It is thus impossible to say 'I want a membrane
# 39 A thick', because if you build such a membrane, it will adapt to its
# inherent thickness during the simulation.  

# Membrane height in A, density should be ~0.861 g/l 
memheight=43.
# Height of the hydrophobic membrane core
memcoreheight=28.
      
# Do we have a target?
if MacroTarget==''
  RaiseError "This macro requires a target. Either edit the macro file or click Options > Macro > Set target to choose a target structure"

Clear
Console off
# Do we already have a scene with water?
waterscene = FileSize (MacroTarget)_water.sce
if waterscene
  LoadSce (MacroTarget)_water
else
  # No scene with water present yet
  # Do we have a scene with the protein embedded in the membrane?
  scene = FileSize (MacroTarget).sce
  if scene
    # Yes, protein with membrane is already present
    LoadSce (MacroTarget)
    # Search for the membrane
    c = CountObj Membrane
    if !c
      RaiseError "The membrane object in (MacroTarget).sce must be named 'Membrane'."
  else
    # No membrane scene present yet
    # Has the user already oriented the protein inside the membrane interactively?
    oriscene = FileSize (MacroTarget)_ori.sce
    if oriscene
      LoadSce (MacroTarget)_ori
      Unselect
    else
      # No orientation scene present yet
      # Load the protein, assuming it's a PDB or YOB file
      for type in 'pdb','yob'
        size = FileSize (MacroTarget).(type)
        if size
          obj = Load(type) (MacroTarget)
          # Align object with major axes to minimize cell size
          NiceOriObj (obj)
          break
      if !size
        RaiseError "Initial structure not found. Make sure to create a project directory and place the structure there in PDB or YOB format"
      Unselect
      # As initial guess, assume the major axis is perpendicular to the membrane
      Style Ribbon
      NiceOri Axis1=Y,Axis2=X
      r = Radius
      AutoMoveTo Z=(r*2),Steps=(delay)
      # Make sure that all hydrogens are present etc.
      Clean
      # Now search for hydrophobic transmembrane helices and strands
      ShowMessage "Searching for transmembrane elements..."
      Wait 1
      # To find exposed hydrophobic residues, we first need the maximum side-chain surface areas
      SurfPar Probe=1.4,Resolution=3,Radii=Solvation
      namelist()='Ala','Val','Leu','Ile','Met','Phe','Tyr','Trp'
      for name in namelist
        obj = BuildRes (name)
        CleanObj (obj)
        AddEnvObj (obj)
        surfmax(name) = SurfAtom Sidechain Obj (obj),Accessible,Unit=Obj
        DelObj (obj)
      # Get a list of hydrophobic side-chain surface areas
      AddEnv
      reslist() = ListRes (join namelist) with contribution to accessible surface of all
      if count reslist
        # Found hydrophobic residues contributing to the surface
        surflist() = SurfAtom Res (join reslist) Sidechain,Accessible,Unit=Res
        # Get a list of all heavily exposed (>30%) hydrophobic side-chains
        for i=1 to count reslist
          name = NameRes (reslist(i))
          if surflist(i)>surfmax(name)*0.3
            # Create list of exposed hydrophobic residues, explist
            explist(count explist+1)=reslist(i)
      # Prepare to calculate the major protein axis, normal to the membrane plane
      axis=0.,0.,0.
      if count explist>5
        # First search for transmembrane helices, then for beta strands (beta barrels)
        for type in 'helix','strand'
          secstrelements=0
          # Get a list of all helices or strands
          atomlist() = ListAtom SecStr (type),compress=yes
          for i=1 to count atomlist step 2
            # Select CA atom of first and last residue
            first = ListAtom CA Res Atom (atomlist(i))
            last = ListAtom CA Res Atom (atomlist(i)+atomlist(i+1)-1)
            totresidues = CountRes Atom (first)-(last)
            hydresidues = CountRes (join namelist) Atom (first)-(last)
            expresidues = CountRes Atom (first)-(last) and (join explist)
            #print '(type) (first)-(last) is (totresidues) residues long, of which (hydresidues) are hydrophobic and (expresidues) are exposed hydrophobic ones' 
            if (type=='helix' and totresidues>17 and hydresidues>7 and expresidues>3) or
               (type=='strand' and totresidues>8 and hydresidues>3 and expresidues>1)
              # Found transmembrane element
              secstrelements=secstrelements+1
              reslist() = ListRes Atom (first) (last),Format='RESName MOLNAME RESNUM'
              ShowMessage 'Found transmembrane (type) (secstrelements) spanning residues (reslist1) to (reslist2)...' 
              MarkAtom (first),(last)
              ColorRes Atom (first)-(last),Yellow
              Wait (delay)
              # Add direction vector to major axis
              px,py,pz,dir() = GroupLine (first)-(last) and CA
              if sum (axis*dir)<0
                # Scalar product is negative, turn direction around
                dir()=-dir
              axis()=axis+dir
              #p1,p2,p3,dir() = GroupLine (first)-(last) and CA,CoordSys=Global
              #ShowArrow Start=Point,(p),End=Point,(p+dir*30)
          if norm axis!=0
            break
        MarkAtom None
        if type=='strand' and secstrelements<5
          ShowMessage "Not enough transmembrane secondary structure elements found that could help to tune the protein orientation."
          axis=0.,0.,0.
          Wait (delay)
      else
        ShowMessage "Not enough exposed hydrophobic residues found. Looking for other features..."
        Wait (delay)
      if not norm axis
        # No transmembrane axis found, look for lipid anchors
        lipid='SmilesMEMBER CH2?CH2CH2CH2CH2?'
        liplist() = ListRes (lipid)
        lipatoms = CountAtom (lipid) 
        for lip in liplist
          id = ListRes (lip),Format="RESNAME RESNUM"
          first,last=SpanAtom Res (lip)
          ShowMessage 'Found transmembrane lipid anchor (id)...' 
          MarkAtom (first),(last)
          ColorAtom (first)-(last),Yellow
          Wait (delay/2)
          px,py,pz,dir() = GroupLine Res (lip) and (lipid)
          if sum (axis*dir)<0
            # Scalar product is negative, turn direction around
            dir()=-dir
          axis()=axis+dir
      if norm axis  
        # Orient the protein such that the axis through the membrane helices is perpenticular to the membrane
        ShowMessage "Orienting transmembrane regions perpendicular to the membrane..."
        oriold() = OriObj 1
        orinew()=oriold
        OriObj 1,0,0,0
        alpha,beta,gamma = OriVec (axis)
        RotateObj 1,Y=(-beta)
        RotateObj 1,Z=(90-gamma)
        # Now rotate around the Y-axis until the user sees most of the protein
        xmax=0
        for i=1 to 30
          x = GroupBox all
          if x>xmax
            xmax=x
            orinew() = OriObj 1
          RotateObj 1,Y=5
        # And rotate to the new orientation
        OriObj 1,(oriold)
        AutoMoveTo X=0,Y=0,Z=(r*3),Steps=(delay),Wait=No
        AutoRotateToObj 1,(orinew),Steps=(delay)
      residuesmax=0
      if count explist or count liplist
        # Scan vertically for highest density of exposed hydrophobic residues or lipid anchors
        ShowMessage "Scanning protein for transmembrane region..."
        plane = ShowPolygon Points,Yellow,Vertices=4, (-r),0,(-r),  (r),0,(-r), (r),0,(r), (-r),0,(r)
        for i=-r to r step 2
          PosObj (plane),Y=(i),Z=(r*3)
          Wait 1
          if count liplist
            # Count lipid atoms in current transmembrane region
            memlipatoms = CountAtom (lipid) GlobalY>(i-memcoreheight/2) GlobalY<(i+memcoreheight/2)
            residues=(0.+memlipatoms)*count liplist/lipatoms
          else
            # Count exposed hydrophobic residues in current transmembrane region
            residues = CountRes (join explist) GlobalY>(i-memcoreheight/2) GlobalY<(i+memcoreheight/2)
          chargedresidues = CountAtom NZ Res Lys or CZ Res Arg or CG Res Asp or CD Res Glu and GlobalY>(i-memcoreheight/2) GlobalY<(i+memcoreheight/2)
          #print 'PosY=(i), ExpResidues=(residues), Charged=(chargedresidues), (residues)'
          residues=residues-chargedresidues*2
          if residues>residuesmax
            residuesmax=residues
            memposy=i
        DelObj (plane)
      # First get the cell that is required to enclose the entire protein, ignoring the membrane
      size1,size2,size3,pos() = GroupBox all,Type=VdW,CoordSys=Global
      if residuesmax
        # Found a membrane embedding
        # Get the cell that is required to enclose the protein
        posmin1()=pos-size*0.5-waterextension
        posmax1()=pos+size*0.5+waterextension
        # Now get the cell that is required to enclose the protein membrane region
        size1,size2,size3,pos() = GroupBox GlobalY>(memposy-memcoreheight/2) GlobalY<(memposy+memcoreheight/2),Type=VdW,CoordSys=Global
        # Along the Y-axis, the size is memheight+waterextension*2
        size2=memheight+waterextension*2-memextension*2
        posmin2()=pos-size*0.5-memextension
        posmax2()=pos+size*0.5+memextension
        # Calculate the new overall cell position, which encloses both cells above
        for i=1 to 3
          for f in 'min','max'
            poslist=pos(f)1(i),pos(f)2(i)
            pos(f)3(i)=(f) poslist
        pos()=(posmin3+posmax3)*0.5
        size()=posmax3-posmin3
        if count explist
          # Color exposed hydrophobic residues orange
          ColorRes (join explist),150
      else
        # No embedding found, put membrane in the middle
        if size2<memheight
          size2=memheight
        size1=size1+memextension*2
        size2=size2+waterextension*2
        size3=size3+memextension*2
        memposy=pos2
        Color Element
        ShowMessage "This protein does not look like a membrane protein. Click 'Continue', then place it in the membrane manually..."
        Wait ContinueButton
      # Create and place cell
      Cell (size)
      PosObj SimCell,(pos)
      if Structure
        # Optimize the hydrogen-bonding network (more stable trajectories)
        ShowMessage "Optimizing hydrogen bonding network..."
        Wait 1
        OptHyd Method=YASARA
      # Create MembPreview, 'MEMBRANE' text that symbolizes the future membrane
      MakeTextObj Name=MembPreview,Width=(size3*10),Height=(size3)
      Font Arial,Height=45%,Color=Grey,Alpha=80,Depth=(memcoreheight),DepthCol=606060,DepthAlpha=50
      PosText 50%,50%,Justify=Center
      Print 'MEMB'
      PosText 50%,0%,Justify=Center
      Print 'RANE'
      ScaleMesh MembPreview,X=((size1/size3)*0.5)
      Style Ribbon
      # Let the membrane preview fly in
      UserInput Off
      ShowMessage "Moving membrane placeholder to suggested position. Exposed hydrophobic residues are colored yellow..."
      PosOriObj MembPreview,Y=-40,Z=-18,Alpha=60
      AutoMoveToObj MembPreview,Y=(memposy),Z=(pos3),Steps=(delay*2),Wait=No
      AutoRotateToObj MembPreview,Alpha=90,Steps=(delay*2)
      UserInput On
      # Let the user tune the embedding
      HUD Obj
      ShowMessage "Move the objects and adjust the cell size to fine-tune the suggested membrane embedding. Click 'Continue' when done."
      Wait ContinueButton
      SaveSce (MacroTarget)_ori
    # Next step is to build the membrane
    # Get the required membrane size
    memsize1,cellsizey,memsize2 = Cell
    if count usermemname
      ShowMessage 'Building a membrane of size (0.0+memsize1) x (0.0+memsize2) A^2 based on user-supplied membrane (usermemname)...'
      membranename='(MacroTarget)_membrane_(0+memsize1)x(0+memsize2)_(usermemname).yob'
    else
      ShowMessage 'Building a membrane of size (0.0+memsize1) x (0.0+memsize2) A^2 with composition (PEApercent)% PEA, (PCHpercent)% PCH, (PSEpercent)%PSE...'
      membranename='(MacroTarget)_membrane_(0+memsize1)x(0+memsize2)_pea(PEApercent)pch(PCHpercent)pse(PSEpercent).yob'
    membrane = FileSize (membranename)
    if membrane
      # Membrane has been built before
      memobj = LoadYOb (membranename)
    else
      # Membrane must be built now
      AutoMoveToObj 1,X=-150,Z=100,Steps=(delay),Wait=No
      AutoMoveToObj MembPreview,X=150,Z=100,Steps=(delay),Wait=No
      AutoMoveToObj SimCell,X=0,Y=150,Z=100,Steps=(delay)
      DelAll
      # Load the equilibrated membrane snapshot that matches best
      if count usermemname
        memobj = LoadYOb (YASARADir)/yob/membrane_(usermemname)
        memfragsize()=usermemsize
      elif PCHpercent>=33
        memobj = LoadYOb (YASARADir)/yob/membrane_pea66pch33
        memfragsize=84.82,82.52
      elif PCHpercent>=20
        memobj = LoadYOb (YASARADir)/yob/membrane_pea80pch20
        memfragsize=85.37,81.02
      else
        memobj = LoadYOb (YASARADir)/yob/membrane_pea
        memfragsize=81.68,79.99
      NameObj 1,Membrane
      PosObj 1,Y=-33,Z=-35
      AutoMoveToObj 1,Z=77,Steps=(delay)
      AutoRotateToObj 1,Alpha=90,Steps=(delay*2),Wait=No
      AutoMoveToObj 1,0,0,(norm memsize*1.5),Steps=(delay*2)
      TransformObj 1
      Cell Auto,Extension=10
      Cell (memfragsize)
      SwitchObj SimCell,off
      # Prepare to mutate membrane residues to get the requested composition.
      # First we search for headgroups that don't have many neighbors and can
      # thus be changed without causing too many bumps. This search is done
      # with periodic boundaries and the original membrane fragment size.
      # Handle both membrane faces separately, each equilibrated membrane fragment consists
      # of 200 residues, residue numbers 1-100 are on side 1, 101-200 on side 2.
      for i=1 to 2
        # Set segment name so that each membrane side can be selected easily afterwards
        siderange='(-99+i*100)-(i*100)'
        SegRes (siderange),MEM(i)
        for name in 'PSE','PCH'
          lipids = CountRes (name) and (siderange)
          if lipids<(name)percent
            # Determine those PEA residues that can best be mutated, i.e. have
            # the smallest number of neighbor atoms that could cause bumps.
            ShowMessage 'Locating exposed lipids on side (i) that can be mutated easily to (name)...'
            Wait 1
            Sim On
            pealist() = ListRes PEA and (siderange),Format=RESNUM
            for j=1 to count pealist
              if name=='PCH'
                neighbors = CountAtom all with distance<5 from N Res (pealist(j))
              else
                neighbors = CountAtom all with distance<5 from PDB 2HA Res (pealist(j)) 
              # Add neighbors*1000 to the residue number so that we can sort easily
              pealist(j)=pealist(j)+neighbors*1000
            Sim Off
            pealist()=sort pealist
            for j=1 to count pealist
              resnum=pealist(j)%1000
              if name=='PCH'
                # Turn phosphatidyl-ethanolamine into phosphatidyl-choline
                SwapAtom H with bond to N Res (resnum),C,Update=Yes
                NameRes (resnum),PCH
              else
                # Turn phosphatidyl-ethanolamine into phosphatidyl-serine
                # SwapAtom creates C35, with an unknown atom number since hydrogens are rebuilt 
                SwapAtom PDB 2HA Res (resnum),C,Update=Yes
                c = ListAtom C35 Res (resnum) 
                NameAtom (c),C
                DelAtom Element H with bond to (c)
                AddHydAtom (c),2
                SwapAtom (c+1) (c+2),O
                Distance (c),all,Bound=Yes,Set=1.231
                SwapBond (c),(c+1) (c+2),1.5
                NameAtom (c+1),OT1
                NameAtom (c+2),OT2
                NameRes (resnum),PSE
              lipids=lipids+1
              ShowMessage 'Created (name) residue (lipids) of ((name)percent) on side (i)...'
              Wait 1
              if lipids==(name)percent
                break
        AutoRotate Y=(180./delay)
        Wait (delay)
        AutoRotate Y=0
      OriObj all,0,0,0
      DelObj SimCell
      # Replicate the membrane fragments 
      axislist='X','Y'
      for i=1 to 2
        memrepsize(i)=memfragsize(i)
        while memrepsize(i)<memsize(i)
          ShowMessage 'Replicating membrane fragment to reach a membrane size of (0+memsize1)x(0+memsize2) A^2.'
          obj = DuplicateObj (memobj)
          MoveObj (obj),(axislist(i))=(memfragsize(i))
          memrepsize(i)=memrepsize(i)+memfragsize(i)
          Wait 1
        JoinObj Membrane,(memobj)
      CenterObj Membrane
      if Atoms>60000
        # Lots of atoms, speed things up
        Antialias off
      AutoMoveToObj Membrane,0,0,Steps=(delay)
      # Keep aspect ratio
      if memrepsize1<memrepsize2
        memrepsize2=memrepsize1*(memsize2/memsize1)
      else
        memrepsize1=memrepsize2*(memsize1/memsize2)
      # Calculate how many lipids we need per side, to 1+ avoids that we round downwards
      lipidsneeded=1+100*(memsize1*memsize2/(memfragsize1*memfragsize2))
      # For each side, iteratively find the region that contains the required number of lipids
      # (We can't do it for both sides together, since we might end up with a different number of lipids per side)
      for i=1 to 2
        print 'Side (i)'
        side='Segment MEM(i)'
        memselmax()=0.5e0*memrepsize
        memselmin()=0.,0.
        lipidsmax = CountRes (side)
        lipidsmin=0
        do
          ShowMessage 'Locating membrane region on side (i) that contains the required number of (lipidsneeded) lipids...'
          Wait 1
          memsel()=(memselmin+memselmax)*0.5
          region='GlobalX>(-memsel1) GlobalX<(memsel1) GlobalY>(-memsel2) GlobalY<(memsel2)'
          lipids = CountRes (side) (region)
          #print 'Min=(memselmin), Max=(memselmax), Sel=(memsel), Lipids=(lipids)'
          if lipids>lipidsneeded
            memselmax()=memsel
            ColorRes (side) and not (region),Grey
          elif lipids<lipidsneeded
            memselmin()=memsel
        while lipids!=lipidsneeded and norm (memselmax-memselmin)>0.1
        # Delete surplus lipids
        DelRes (side) and not (region)
      RenumberRes all,1  
      # Enclose in cell, which is now larger than needed
      Cell Auto,Extension=0
      _,_,z = Cell
      Cell Z=(z+20)
      # Shrink the cell step-by-step while minimizing the lipids
      Longrange None
      Cutoff 5.24
      Temp 298K
      # Compressed membrane size is 5% smaller to close holes faster
      f=0.95
      commemsize()=memsize*f
      TimeStep 2,1
      Sim On
      # Since this is done in vacuo, prevent electrostatic interactions
      ChargeAtom Obj Membrane,0
      shrinkstep=0.6e0
      while 1
        x,y = Cell
        ShowMessage 'Compressing membrane to reach a size of (0+commemsize1)x(0+commemsize2) A^2, current size is (0+x)x(0+y) A^2...'
        if x<=commemsize1 and y<=commemsize2
          break
        Sim Off
        if x>commemsize1
          x=x-shrinkstep
        if y>commemsize2
          y=y-shrinkstep
        Cell X=(x),Y=(y)
        Sim On
        TempCtrl SteepDes
        Wait SpeedMax<3000
        TempCtrl Rescale
        for i=1 to 10
          StabilizeMembrane 'Z',35.,0
          Wait 10,femtoseconds
      # Allow membrane to relax
      ShowMessage 'Running short MD simulation of the compressed membrane to close gaps...'
      Sim On
      TempCtrl Anneal
      Wait 200,femtoseconds
      TempCtrl Rescale
      Temp 298K
      for i=1 to 100
        StabilizeMembrane 'Z',35.,0
        Wait 10,femtoseconds
      TempCtrl Anneal
      Wait 200,femtoseconds
      # Uncompress the membrane
      ScaleObj Membrane,X=(1./f),Y=(1./f)
      # The membrane ends have now been fused, the fusion region is at the periodic boundary.
      # Since the fusion region is not equilibrated and less ideally packed, we shift the
      # membrane by half the cell size to make sure that the worst part ends up in the cell
      # center - where it will be deleted anyway, since it overlaps with the protein.
      MoveAtom all,(memsize*0.5)
      Cell (memsize)
      # Wrap the lipids
      Sim On
      Sim Off
      # Finished
      Color Element
      SaveYOb 1,(membranename)
      AutoMoveToObj 1,Z=-45,Steps=(delay)
      # Load the orientation scene again and include the membrane
      LoadSce (MacroTarget)_ori
      memobj = LoadYOb (membranename)
    # Replace the 'MEMBRANE' text with the real membrane
    NameObj (memobj),Membrane
    TransferObj Membrane,MembPreview,Local=Keep
    SwapObj MembPreview,Membrane
    DelObj MembPreview
    # Align the cell with the global coordinate system in case it has been rotated by the user
    alpha,beta,gamma = OriObj SimCell
    Rotate Y=(-beta)
    Rotate Z=(-gamma)
    Rotate X=(-alpha)
    # The membrane placeholder and the actual membrane have top and bottom flipped
    RotateObj Membrane,X=180
    # Use the opportunity to remember the local Y-coordinate of the membrane center in the cell
    _,memposy = PosObj Membrane
    _,cellposy = PosObj SimCell
    memposy=memposy-cellposy+cellsizey*0.5
    # Freeze the current orientation
    TransformObj all
    # Make sure that protein and membrane share the same local coordinate system
    # The protein origin 0/0/0 is then at the center of the membrane
    TransferObj 1,Membrane,Local=Fix
    # Calculate the bounding box of the transmembrane region
    transmemsize1,_,transmemsize2 = GroupBox Obj 1 LocalY>(-memcoreheight/2+2) LocalY<(memcoreheight/2-2), Type=Nuclear
    #ColorAtom Obj 1 LocalY>(-memcoreheight/2+2) LocalY<(memcoreheight/2-2),magenta
    # The following minimizations will be done quickly
    Cutoff 5.24
    Longrange None
    # Calculate initial XZ-scaling factor required to shrink the membrane region
    # 12 was too much, 6 was too little
    scalexz=(max transmemsize-10)/max transmemsize
    # Prepare to embed the protein in the membrane. We cannot simply delete bumping
    # lipids: if a single overlap was enough to delete a lipid, we would end up with too few 
    # lipids and large holes between membrane and protein. Instead we first shrink the protein,
    # then delete bumping lipids, and finally grow the protein again during an energy minimization,
    # so that it can slowly push the surrounding lipids away until they smoothly cover the protein.
    #
    # Initialize simulation before the object is scaled, otherwise the incorrect
    # bond lengths prevent YASARA from caching the parameters to be on the safe side
    Sim Init
    # First shrink the protein strongly, in case it has a large pore
    ScaleObj 1,X=(scalexz*0.5),Z=(scalexz*0.5)
    ShowMessage "Deleting lipids that strongly bump into the shrunk protein..."
    Wait 1
    # Delete bumping lipids
    DelRes Obj Membrane with distance<0.75 from Obj 1
    ScaleObj 1,X=2,Z=2
    DelRes Obj Membrane with distance<0.75 from Obj 1
    # During the expansion, we don't want the protein atoms to move
    FixObj 1
    scalexzstart=scalexz
    ShowMessage "Expanding protein to fill the membrane pore..."
    Sim On
    ChargeAtom Obj Membrane,0
    # Iterate while the protein has not reached its original size yet
    while scalexz<1
      # Quick energy minimization
      Sim On
      TempCtrl SteepDes
      Wait SpeedMax<4000
      Temp 298K
      TempCtrl Rescale
      for i=1 to 20
        StabilizeMembrane 'Y',memposy,0
        Wait 10,femtoseconds
      f=1.02
      ScaleObj 1,X=(f),Z=(f)
      ShowMessage 'Expanding protein to fill the membrane pore, (0+(scalexz-scalexzstart)*100/(1.-scalexzstart))% completed...'
      scalexz=scalexz*f
    # Set original size again, since we did not reach 1.0 exactly
    ScaleObj 1,X=(1./scalexz),Z=(1./scalexz)
    ShowMessage "Short energy minimization to optimize membrane geometry..."
    # Minimize before saving, since the pulling caused incorrect bond lengths
    Sim On
    TempCtrl Anneal
    Wait 400,femtoseconds
    Sim Off
    SaveSce (MacroTarget)
  # Set reliable simulation parameters again
  Cutoff (cutoff)
  Longrange (longrange)
  # Get membrane Y-position in the local coordinate system of the simulation cell
  # (needed in case we have not just created the membrane but loaded it instead).
  Sim On
  _,memposy = PosAtom Obj Membrane,Mean=yes
  Sim Off
  # Fill the cell with water
  Free
  Experiment Neutralization
    WaterDensity 0.998
    pH (ph)
    NaCl (nacl)
    pKaFile (MacroTarget).pka
    #Vacuum LocalY>(memposy-memcoreheight*0.5) LocalY<(memposy+memcoreheight*0.5)
    Speed Fast
  Experiment On
  Wait ExpEnd
  # Delete leftover waters in the membrane region, except those that are far away
  # from the membrane (waters inside the channel of a pore protein) 
  ShowMessage 'Deleting water molecules inside the membrane...'
  Wait 1
  Sim On
  memregion='LocalY>(memposy-memcoreheight*0.5) LocalY<(memposy+memcoreheight*0.5)'
  DelRes HOH Atom O (memregion) with distance<9 from Obj Membrane
  # Tag those waters that are allowed in the membrane region (usually pore waters)
  # No pulling force will be applied to these waters during the equilibration period
  Sim On
  PropRes all,0
  PropRes HOH Atom O (memregion),100
  Sim Off
  StickRes Water
  # Save scene with water
  SaveSce (MacroTarget)_water
  HideMessage
Wait 1

if constrain=='yes'
  # Constrain bond lengths to hydrogens
  FixBond all,Element H
  # Constrain bond angles in water
  FixAngle Water,Water,Water
  # Multiple timestep: 1.3333 femtoseconds for intramolecular and 3*1.3333 = 4 fs for intermolecular forces
  timestep=3,1.3333
  # Save simulation snapshots every 6250 simulation steps
  # (with a timestep of 4 femtoseconds, that's 6540*4 fs = 25 picoseconds).
  savesteps=6250
else
  # No constraints
  FreeBond all,all
  FreeAngle all,all,all
  # Smaller timestep, since we don't use constraints: 2*1.25 = 2.5 fs
  timestep=2,1.25
  # Save simulation snapshots every 10000 simulation steps
  # (with a timestep of 2.5 femtoseconds, that's 10000*2.5 fs = 25 picoseconds).
  savesteps=10000

# Set final simulation parameters
Temp (temperature)
Cutoff (cutoff)
Longrange (longrange)
TimeStep (timestep)
# Make sure all atoms are free to move
Free
# Alread a snapshot/trajectory present?
i=00000
filename='(MacroTarget)(i).(format)'
running = FileSize (filename)
if !running
  # First an energy minimization
  Experiment Minimization
  Experiment On
  Wait ExpEnd
  # And now start with the real simulation
  Sim On
else
  # Simulation has been running before
  ShowMessage "Simulation has been running before, loading last snapshot..."
  Wait 1
  # Switch console off to load the snapshots quickly
  if format=='sim'
    # Find and load the last 'sim' snapshot
    do
      i = i+1
      found = FileSize (MacroTarget)(i).sim
    while found
    LoadSim (MacroTarget)(i-1)
  else
    # XTC format requires that the entire trajectory is read in to find the last one
    do
      i = i+1
      eof,time = LoadXTC (filename),(i)
      ShowMessage 'Searching XTC trajectory for last snapshot, showing snapshot (i) at (0+time) fs'
      Sim Pause
      Wait 1
    while !eof
    Sim Continue
HideMessage

# Set pressure control, we combine a manometer with solvent density measurements
PressureCtrl Combined,Pressure=(pressure),Name=(solvent),Density=(density)

# And finally, make sure that future snapshots are saved
Save(format) (filename),(savesteps)

# At the beginning of the simulation, the membrane is not ideally packed yet and needs
# about 250ps equilibration time. During this period it is still 'vulnerable' to
# water molecules that are squeezed in. We thereforce keep these waters out.
t = Time
if t<equiperiod*1000
  # The membrane stabilization forces add kinetic energy, strengthen temperature coupling
  TempCtrl Rescale,RelaxTime=100
  ShowMessage 'Running (equiperiod)ps of equilibration simulation, preventing water from entering the membrane...'
  while t<equiperiod*1000
    # Get current membrane Y position
    _,memposy = PosAtom Obj Membrane,Mean=yes
    for i=1 to 10
      StabilizeMembrane 'Y',memposy,'Water'
      Wait 20,femtoseconds
    t = Time
  HideMessage
# Set final temperature control
TempCtrl Rescale,RelaxTime=100

if duration=='forever'
  Console On
  Wait forever
else
  # Wait for given number of picoseconds
  do
    Wait 10
    t = Time
  while t<duration*1000
  Sim Off

# STABILIZE A MEMBRANE BY PULLING THE LIPIDS 
# ==========================================
# 'axis' is the membrane axis (normally 'Y' or 'Z'), 'pos' is the
# position of the membrane center along 'axis'. 'water' is either 0
# (no water) or selects the water object to be kept out of the membrane.
def StabilizeMembrane axis,pos,water
  global memheight,memcoreheight
  
  # Set single atom force f and residue acceleration a
  s = SimSteps
  f=3000.*2*s
  a=10000.*2*s
  # Accelerate complete lipid residues that try to leave the membrane
  AccelRes Obj Membrane Local(axis)>(pos+memheight*0.5),(axis)=(-a)
  AccelRes Obj Membrane Local(axis)<(pos-memheight*0.5),(axis)=(a)
  # Pull up the lipid heads that get sucked into the membrane
  # Replace 'P' with the name of the phosphorus atom in your lipid if needed
  ForceAtom P Segment MEM1 Local(axis)>(pos-14),(axis)=(-f)
  ForceAtom P Segment MEM2 Local(axis)<(pos+14),(axis)=(f)
  # Pull back the lipid tails that get too close to the membrane surface
  # Replace 'C34' and 'C18' with the names of the terminal carbons in the long
  # and short lipid tails
  ForceAtom C34 Segment MEM1 Local(axis)<(pos+2),(axis)=(f)
  ForceAtom C18 Segment MEM1 Local(axis)<(pos-2),(axis)=(f)
  ForceAtom C34 Segment MEM2 Local(axis)>(pos-2),(axis)=-(f)
  ForceAtom C18 Segment MEM2 Local(axis)>(pos+2),(axis)=-(f)
  if water
    # Keep the water out of the membrane (atoms with property 100 are allowed membrane pore waters)
    AccelRes Property!=100 Local(axis)>(pos-memcoreheight*0.5) Local(axis)<(pos) Obj (water),(axis)=(-a)
    AccelRes Property!=100 Local(axis)<(pos+memcoreheight*0.5) Local(axis)>(pos) Obj (water),(axis)=(a)
    # Ions need more force
    AccelRes CIP CIM Local(axis)>(pos-memcoreheight*0.5) Local(axis)<(pos) Obj (water),(axis)=(-a*10)
    AccelRes CIP CIM Local(axis)<(pos+memcoreheight*0.5) Local(axis)>(pos) Obj (water),(axis)=(a*10)
    
    Color Element
    ColorRes Property!=100 Local(axis)>(pos-memcoreheight*0.5) Local(axis)<(pos) Obj (water),Magenta
    ColorRes Property!=100 Local(axis)<(pos+memcoreheight*0.5) Local(axis)>(pos) Obj (water),Grey
  
  """
  # Coloring
  Color Element
  ColorRes Obj Membrane Local(axis)>(pos+memheight*0.5),Green
  ColorRes Obj Membrane Local(axis)<(pos-memheight*0.5),Yellow
  ColorAtom P Segment MEM1 Local(axis)>(pos-14),Red
  ColorAtom P Segment MEM2 Local(axis)<(pos+14),Red
  ColorAtom C34 Segment MEM1 Local(axis)<(pos+2),Magenta
  ColorAtom C18 Segment MEM1 Local(axis)<(pos-2),Magenta
  ColorAtom C34 Segment MEM2 Local(axis)>(pos-2),Blue
  ColorAtom C18 Segment MEM2 Local(axis)>(pos+2),Blue
  """