#! /usr/bin/python3 -t 

MODULE_NAME = "mformula_oxocarbon"
MODULE_DESC = "general routines for building oxocarbon formulas"
MODULE_VERS = "1.0"

MODULE_COPYRIGHT = "Copyright © 2009-07-11 by the State University of Campinas (UNICAMP)"

import sys
import os
import re
import copy
import math; from math import sqrt,sin,cos,pi
# sys.path[1:0] = [ sys.path[0] + '/../lib', os.path.expandvars('${STOLFIHOME}/lib'), '.' ]
# sys.stderr.write(re.sub('[,]', ',\n', "%s: path = %r\n" % (PROG_NAME, sys.path)))
import rn
import mformula

def build_formula_C_n_O(svg,n,v0,v1,chL,chR,chD) :
  """Builds the formula for a linear C[n]O molecule.
  
  The molecule is linear, with the O atom at right. Even-numbered bonds will have valency {v0}, odd-numbered bonds will have valency {v1}. The first carbon will have a charge {chL/chD}. The oxygen will have a charge {chR/chD}."""
  fm = mformula.obj()
  
  xa = 0         # Position of first atom. 
  ka = fm.add_atom("C", (xa, 0), chL, chD)
  
  xk = xa  # {X} coord of previous atom.
  vk = v0  # Valence of bond to next atom.
  for k in range(1,n) :
    # Add atom number {k} (a carbon):
    xk += svg.rel_bond_length(vk)
    kk = fm.add_atom("C", (xk, 0), 0,1)
    fm.add_bond(kk-1,kk,vk)
    # Compute valency of next bond:
    vk = [ v0, v1 ] [k % 2]
    
  # Add the final oxygen:
  xO = svg.CObondlength * svg.rel_bond_length(vk)
  kO = fm.add_atom("O", (xO, 0), chR,chD)
  assert kO == n

  # Bind it:
  fm.add_bond(kO-1,kO,vk)

  if (svg.style[0] == 'S') and (chL != 0) :
     # Add charge symbol on first carbon:
     paq = (xa, svg.CHbondlength)
     kaq = fm.add_atom("+", paq, chL,chD)
 
  if (svg.style[0] == 'S') and (chR != 0) :
     # Add charge symbol on oxygen:
     pOq = (xO, svg.CHbondlength)
     kOq = fm.add_atom("+", pOq, chR,chD)

  return fm
  #----------------------------------------------------------------------
  
def build_formula_C_n_O_2(svg,n,v0,v1,chL,chR,chD) :
  """Builds the formula for a linear OC[n]O molecule.
  
  The molecule is linear, with the O atom at right. Even-numbered bonds
  will have valency {v0}, odd-numbered bonds will have valency {v1}. The
  first oxygen will have a charge {chL/chD}. The last oxygen will have a
  charge {chR/chD}."""
   
  fm = mformula.obj()
    
  xO1 = 0
  # Add the initial oxygen:
  kO1 = fm.add_atom("O", (xO1,0),chL,chD)
  assert kO1 == 0
  
  xk = xO1            # Coordinate of previous atom.
  vk = v0             # Valency of bond to next atom.
  for k in range(n) :
    # Compute the position of the next atom:
    xk = xk + (svg.CObondlength if k == 0 else 1.0)*svg.rel_bond_length(vk)
    # Add atom number {k+1} (a carbon):
    kk = fm.add_atom("C", (xk,0), 0,0)
    assert kk == k + 1
    # Bind it:
    fm.add_bond(kk-1,kk,vk)
    # Compute valency of next bond:
    vk = [ v0, v1 ] [(k+1) % 2]
    
  # Compute the position of the final oxygen:
  xO2 = xk + svg.CObondlength*svg.rel_bond_length(vk)
  # Add the final oxygen:
  kO2 = fm.add_atom("O", (xO2,0),chR,chD)
  assert kO2 == n+1

  # Bind it:
  fm.add_bond(kO2-1,kO2,vk)

  if (svg.style[0] == 'S') and (chL != 0) :
     # Add charge symbol on left oxygen:
     p0q = (xO1, svg.CHbondlength)
     k0q = fm.add_atom("+", p0q, chL,chD)

  if (svg.style[0] == 'S') and (chR != 0)  :
     # Add charge symbol on right oxygen:
     pnq = (xO2, svg.CHbondlength)
     knq = fm.add_atom("+", pnq, chR,chD)

  return fm
  #----------------------------------------------------------------------
  
def build_formula_cyclic_CO_n(svg,n,ch,subs) :
  """Builds the formula for a cyclic polyketone (OC)[n] molecule or anion.
  
  The molecule is cyclic; the bottom side of the ring is horizontal. The
  parameter {ch} is the total charge and should be an integer in
  {-n..0}. The valencies are computed from {n} and {ch}. All bonds have
  the same length which depends on the mean valency around the ring. The
  parameter {subs} is a string with {n} letters, describing the
  'oxygen-like' groups CCW around the ring; each letter is 'O' for
  oxygen, 'D' for dicyanomethylene. The carbons corresponding to the
  first and last letters of {subs} define a horizontal bond at the
  bottom of the ring."""

  fm = mformula.obj()
  assert n >= 2
  assert n == len(subs)
  
  kk = [None] * n # {kk[i]} is the C atom of group {i}.
  
  # Compute the mean charge {q} on the oxygen atoms:
  q = (ch + 0.0) / n
  
  # Compute the valency {ccv} and length {ccb} of the C-C ring bond:
  ccv = 1.0 - q/2.00 
  if n == 2 :
    # Ring bonds become doubled:
    ccb = svg.rel_bond_length(2*ccv)
  else :
    ccb = svg.rel_bond_length(ccv)  # Length of bonds along main ring.
  sys.stderr.write("! build_formula_cyclic_CO_n  n = %r  ch = %r  q = %r  ccv = %r  ccb = %r\n" % (n,ch,q,ccv,ccb))
     
  R = (ccb/2.00)/sin(pi/n) # Radius of ring (mol. center to atom center).

  # Build a formula {fg} with one carbonyl group at top:
  fgO = build_formula_cyclic_CO_elem(svg,R, ch,n)
  fgD = build_formula_cyclic_CCCN2_elem(svg,R, ch,n)

  # Add the carbonyl groups:
  for i in range(n) :
    sys.stderr.write("i = %d  subs[i] = %r\n" % (i,subs[i]))
    if subs[i] == 'O' :
      fg = fgO
    elif subs[i] == 'D' :
      fg = fgD
    else :
      assert False
    kk[i] = fm.add_subformula(fg, 0, ((i+0.5)/n + 0.5)*360.0, [0,0])
  
  if n == 2 :
    # Add one bond with double the valency {ccv}:
    k0 = kk[0]
    k1 = kk[1]
    fm.add_bond(k0+0,k1+0,2*ccv)
  elif n > 2 :
    # Add bonds around the ring with valency {ccv}:
    for i in range(n) :
      k0 = kk[i]
      k1 = kk[(i+1) % n]
      fm.add_bond(k0+0,k1+0,ccv)

  return fm
  #----------------------------------------------------------------------
  
def build_formula_cyclic_CO_n_acid(svg,n,subs,ROT) :
  """Builds the formula for a cyclic polyketo diol (OC)[n-2](COH)[2] molecule.
  
  The molecule is cyclic. The parameter {subs} is a string with {n}
  letters, describing the 'oxygen-like' groups CCW around the ring; each
  letter is 'H' or 'h' for hydroxyl, 'O' for oxygen, 'D' for
  dicyanomethylene. There must be an equal number of 'H's and 'h's and
  they must occur as consecutive pairs 'Hh'. Carbons which bear an 'H'
  are connected by a double bond to the following 'h'-bearing carbon.
  All bonds have the same length which depends on the mean valency. The
  whole molecule is rotated by {ROT} degrees CCW (but the OHs are
  horizontal in 'S' style); if {ROT} is zero, the carbons corresponding
  to the first and last letters of {subs} define a horizontal bond at
  the bottom of the ring."""

  assert n >= 3
  assert n == len(subs)

  # Build the full formula {fm}:
  fm = mformula.obj()
  
  kk = [None] * n # {kk[i]} is the C atom of group {i}.
  dk = [None] * n # {dk[i]} is 1 iff the bound to the previous atom is double.
  
  # Count the number of 'H's in the string:
  nH = 0; nh = 0
  for i in range(len(subs)) :
    if subs[i] == 'H' :  nH = nH + 1
    if subs[i] == 'h' :  nh = nh + 1
  assert nH == nh

  # Compute the valency {ccv} and length {ccb} of the C-C ring bond:
  ccv = (n + nH + 0.0)/(n + 0.0) 
  ccb = svg.rel_bond_length(ccv)
  sys.stderr.write("! build_formula_cyclic_CO_n_acid  n = %r  nH = %r  nh = %r  ccv = %r  ccb = %r\n" % (n,nH,nh,ccv,ccb))
     
  R = (ccb/2.00)/sin(pi/n) # Radius of ring (mol. center to atom center).

  # Build common subformulas:
  fgO = build_formula_cyclic_CO_elem(svg,R, 0,1)
  fgD = build_formula_cyclic_CCCN2_elem(svg,R, 0,1)

  # Add the carbons + side groups:
  for i in range(n) :
    sys.stderr.write("i = %d  subs[i] = %r\n" % (i,subs[i]))
    # Compute angle {ae} of element axis (degrees, CCW from 12:00):
    ae = ((i+0.5)/n - 0.5)*360.0 + ROT
    # Reduce to range -180 to +180:
    while ae < -180 : ae = ae + 360
    while ae > +180 : ae = ae - 360
    # Pick the right element {fg}:
    if subs[i] == 'O' :
      fg = fgO
    elif subs[i] == 'D' :
      fg = fgD
    elif subs[i] == 'H' or subs[i] == 'h' :
      # Compute bend angle {ah} of OH bond (degrees, CCW from elem axis) so that the OH is horizontal:
      if svg.style[0] == 'B' :
        ah = 0
      elif svg.style[0] == 'S' :
        if ae < +0.1 : 
          ah = -ae - 90
        elif ae > +179.9 : 
          ah = 270 - ae
        else :
          ah = 90 - ae
      else :
        assert False
      fg = build_formula_cyclic_COH_elem(svg,R,ah)
    else :
      assert False
    # sys.stderr.write("! build_formula_cyclic_CO_n_acid  i = %d  xflip = %r\n" % (i, xflip))
    kk[i] = fm.add_subformula(fg, 0, ae, [0,0])
  
  # Decide valence of ring bonds: {cch} between 'H' and 'h' and {ccx} of other pairs:
  if nH + nh == n :
    cch = 1.5; ccx = 1.5
  else :
    cch = 2.0; ccx = 1.0
  
  # Add bonds around the ring with valency {ccv}:
  for i0 in range(n) :
    i1 = (i0+1) % n
    k0 = kk[i0]
    k1 = kk[i1]
    # sys.stderr.write("! build_formula_cyclic_CO_n_acid  subs[%d] = %r  subs[%d] = %r\n" % (i0,subs[i0],i1,subs[i1]))
    if subs[i1] == 'h' :
      assert subs[i0] == 'H'
      fm.add_bond(k0+0,k1+0,cch)
    else :
      assert subs[i0] != 'H' 
      fm.add_bond(k0+0,k1+0,ccx)

  return fm
  #----------------------------------------------------------------------
  

def build_formula_cyclic_CO_elem(svg,R,chN,chD) :
  """Builds the formula for one CO unit of a cyclic polyketone molecule or anion.
  
  The group is vertical with the carbon (atom 0) centered at {[0,R]}.
  The oxygen will have a charge {chN/chD}. The C-O valency is computed
  from {chN/chD}."""
   
  fm = mformula.obj()
  # Compute the valency {vco} of the C-O bond:
  vco = 2.0 + chN/chD
  
  yC = R
  iC = fm.add_atom("C", (0, yC), 0,0)
  
  yO = yC + svg.CObondlength*svg.rel_bond_length(vco)
  iO = fm.add_atom("O", (0, yO), chN,chD)

  if (svg.style[0] == 'S') and (chN != 0) :
     # Add charge symbol on oxygen:
     yq = yO + svg.OHbondlength
     kq = fm.add_atom("+", (0, yq), chN,chD)
  
  fm.add_bond(iC,iO,vco)

  return fm
  #----------------------------------------------------------------------
  
def build_formula_cyclic_CCCN2_elem(svg,R,chN,chD) :
  """Builds the formula for one CC(CN)2 unit of a cyclic pseudo-polyketone 
  molecule or anion.
  
  The group is vertical with the carbon (atom 0) centered at {[0,R]}.
  The central carbon will have a charge {chN/chD}. The valency of the
  root C-C bond is computed from {chN/chD}."""
   
  fm = mformula.obj()

  # Add atoms:
  pC = (0, R)
  iC = fm.add_atom("C", pC, 0,0) # Root carbon.
  
  q = chN/chD

  vCM = 2.0 + q
  dCM = svg.rel_bond_length(vCM)
  pM = (0, pC[1] + dCM)
  iM = fm.add_atom("C", pM, chN,chD) # Middle carbon.

  # Distance between carbon M and cyanide carbons:
  dCS = svg.rel_bond_length(1.0)
  
  # Create a cyanide group:
  fg = build_formula_cyanide_elem(svg,dCS)
  
  # Choose the tilt angle of the cyanide groups:
  a0 = 60*pi/180
  a1 = pi - a0
  iR = fm.add_subformula(fg, 0, a0*180/pi, pM)
  iS = fm.add_subformula(fg, 0, a1*180/pi, pM)

  if (svg.style[0] == 'S') and (chN != 0) :
     # Add charge symbol:
     pMq = (pM[0]+svg.CHbondlength, pM[1])
     kMq = fm.add_atom("+", pMq, chN,chD)

  # Add bonds:
  fm.add_bond(iM,iR,1.0)
  fm.add_bond(iM,iS,1.0)
  fm.add_bond(iC,iM,vCM)
  return fm
  #----------------------------------------------------------------------
  
def build_formula_cyanide_elem(svg,R):
  """Build a cyanide group (horizontal, carbon abscissa {R} of origin):"""
  fg = mformula.obj()
  tb = svg.rel_bond_length(3.0)
  iA = fg.add_atom("C", (R, 0), 0,0)
  iB = fg.add_atom("N", (R + tb, 0),  0,0)
  fg.add_bond(iA,iB,3.0)
  return fg
  #----------------------------------------------------------------------

def build_formula_cyclic_COH_elem(svg,R,ang) :
  """Builds the formula for one COH unit of a cyclic pseudo-polyketone
  acid molecule or anion.
  
  The group is vertical with the carbon (atom 0) centered at {[0,R]} and
  the hydrogen tilted by {ang} degrees counterclockwise."""
   
  fm = mformula.obj()

  vCO = 1.0                 # Valence of C-O bond.
  dCO = svg.CObondlength*svg.rel_bond_length(vCO)  # Length of C-O bond.
  
  vOH = 1.0              # Valence of O-H bond.
  dOH = svg.OHbondlength  # Length of O-H bond.
  
  # Choose the C-O-H bend angle:
  a0 = ang*pi/180
  ca0 = cos(a0)
  sa0 = sin(a0)

  # Points:
  pC = [0, R]                           # Root carbon.
  pO = rn.add(pC, [0, dCO])             # Oxygen.
  pH = rn.add(pO, [-sa0*dOH, +ca0*dOH]) # Hydrogen.
  
  # Add atoms:
  iC = fm.add_atom("C", pC, 0,0) # Root carbon.
  iO = fm.add_atom("O", pO, 0,0) # Oxygen.
  iH = fm.add_atom("H", pH, 0,0) # Hydrogen.
  
  # Add bonds:
  fm.add_bond(iC,iO,vCO)
  fm.add_bond(iO,iH,vOH)
  return fm
  #----------------------------------------------------------------------

def build_formula_3_4_dialkynyl_3_cyclobutene_1_2_dione(svg,n,aleg) :
  """Builds the formula for a 3_4_dialkynyl_3_cyclobutene_1_2_dione biradical.
  
  The radical is shaped for assembling a cyle with {n} units centered at
  [0,0]. The bonds exit the cyclobutene ring at an angle {aleg} with the
  double bond. The biradical is symmtric around the Y axis, with the
  cyclobutene at left."""

  fm = mformula.obj()

  dCC1 = svg.rel_bond_length(1.0) # Relative length of single C-C bonds.
  dCC2 = svg.rel_bond_length(2.0) # Relative length of double C-C bonds.
  dCO2 = svg.CObondlength*svg.rel_bond_length(2.0) # Relative length of double C-O bonds.
  tb = svg.rel_bond_length(3.0) # Relative length of triple bonds.

  abend = (pi/n - aleg)/2.0  # External bend angle at the acetenyl carbons.
  
  # The cyclobutene angles are distorted by the double bond:
  a0 = math.asin(0.5 - dCC2/2.0)  # Lift angle from top right to top left cyclobutene carbons.
  sys.stderr.write("a0 = %.2f\n" % (a0*180/pi))
  ca0 = cos(a0)
  sa0 = sin(a0)
  
  # The cyclobutene-acetenyl links exit as specified:
  a1 = pi/2 - aleg       # Lift angle from top left cyclobutene carbon to first acetenyl carbon.
  sys.stderr.write("a1 = %.2f\n" % (a1*180/pi))
  ca1 = cos(a1)
  sa1 = sin(a1)
  
  # Choose direction of C=O in carbonyls to bisect the corner angle:
  a2 = (a0 + pi/2)/2.0     # Lift angle of C=O bond.
  sys.stderr.write("a2 = %.2f\n" % (a2*180/pi))
  ca2 = cos(a2)
  sa2 = sin(a2)

  # The first single-triple bonds are bent as needed: 
  a3 = a1 - abend             # Lift angle of triple C-C bond in acetenyl.
  sys.stderr.write("a3 = %.2f\n" % (a3*180/pi))
  ca3 = cos(a3)
  sa3 = sin(a3)

  # The mid-leg single-triple bonds are bent as needed, too: 
  a4 = a3 - abend             # Lift angle of single bond between acetenyls.
  sys.stderr.write("a3 = %.2f\n" % (a3*180/pi))
  ca4 = cos(a4)
  sa4 = sin(a4)
  
  # Atom positions, with origin at center of double bond:
  p0 = [ +dCC2*0.0/2, +dCC2*1.0/2 ]         # Top right carbon of cyclobutene.
  p1 = rn.add(p0, [ -dCC1*ca0, +dCC1*sa0 ]) # Top left carbon of cyclobutene.
  p2 = rn.add(p1, [ -dCO2*ca2, +dCO2*sa2 ]) # Oxygen of top left carbonyl.
  p3 = rn.add(p0, [ +dCC1*ca1, +dCC1*sa1 ]) # Bottom carbon of acetenyl.
  p4 = rn.add(p3, [ +tb*ca3, +tb*sa3 ]) # Top carbon of acetenyl.

  # Displacement {ux,uy} from middle of double bond to join between biradicals:
  ux = p4[0] + dCC1*ca4/2.0
  uy = p4[1] + dCC1*sa4/2.0
  
  # Displacement to bring the molecular center to the origin:
  dsp = [ -(ux + uy/math.tan(pi/n)), 0.0 ]
  
  # Build a formula {fh} for the top half of the biradical:
  fh = mformula.obj()
  k0 = fh.add_atom("C", rn.add(p0, dsp), 0,0)
  k1 = fh.add_atom("C", rn.add(p1, dsp), 0,0)
  k2 = fh.add_atom("O", rn.add(p2, dsp), 0,0)
  k3 = fh.add_atom("C", rn.add(p3, dsp), 0,0)
  k4 = fh.add_atom("C", rn.add(p4, dsp), 0,0)
  fh.add_bond(k0,k1,1.0)
  fh.add_bond(k1,k2,2.0)
  fh.add_bond(k0,k3,1.0)
  fh.add_bond(k3,k4,3.0)

  # Build the full biradical with two copies of the half one:
  m0 = fm.add_subformula(fh, 0,   0,  [0,0])
  m1 = fm.add_subformula(fh, 1, 180,  [0,0])
  fm.add_bond(m0+0,m1+0,2.0)
  fm.add_bond(m0+1,m1+1,1.0)
  return fm
  #----------------------------------------------------------------------
 
def build_formula_COO_1neg(svg,ang) :
  """Builds the formula for a carboxylate group, with distributed -1 charge.
  
  The carbon is atom number 0 and is at the origin. The carboxylate is
  Y-symmetrical and points to the right. The parameter {ang} is the
  angle (in degrees) between the oxygens. """

  fm = mformula.obj()
  
  # Compute angle {t} (degrees) of oxygens from vertical:
  t = (180 - ang)/2.0
  ct = cos(t*pi/180)
  st = sin(t*pi/180)
  
  # Bond valencies and lengths:
  oq = -1.0/2.0  # Mean charge of O atom.
  ov = 2.0 + oq  # Mean valence of C-O bond.
  ob = svg.rel_bond_length(ov) # Relative length of C-O bond.

  # Atom centers:
  p0 = [0,0]                            # Carbon.
  p1 = rn.add(p0, rn.scale(ob, [+st, -ct]))   # Lower oxygen.
  p2 = rn.add(p0, rn.scale(ob, [+st, +ct]))   # Upper oxygen.

  # Place atoms:
  k0 = fm.add_atom("C", p0,  0,0)
  k1 = fm.add_atom("O", p1, -1,2)
  k2 = fm.add_atom("O", p2, -1,2)

  if svg.style[0] == 'S' :
     # Add charge symbol:
     pq = rn.add(p0, rn.scale(ob, [0.65 + st, 0]))
     kq = fm.add_atom("+", pq, -1,1)

  # Insert bonds (with dashes on distal side):
  fm.add_bond(k0,k1,ov)
  fm.add_bond(k2,k0,ov)
  
  return fm
  #----------------------------------------------------------------------
 
def build_formula_COOH(svg,AOO,AOH,eq) :
  """Builds the formula for a carboxyl group -(CO)OH.
  
  The carbon is atom number 0 and is at the origin. The carboxylate
  points to the right. The parameter {AOO} is the angle (in degrees)
  between the oxygens. If {eq} is zero the hydrogen is attached to the
  lower oxygen, with a bend angle {AOH} (degrees, CCW). If {eq} is 1 the
  hydrogen is sort of halfway between the two oxygens and {AOH} is
  ignored."""

  fm = mformula.obj()
  
  # Compute dip angle {ao1} (radians) of the lower oxygen:
  ao1 = 0.5*AOO*pi/180
  cao1 = cos(ao1)
  sao1 = sin(ao1)
  
  # Compute rise angle {ao2} (radians) of the upper oxygen:
  ao2 = AOO*pi/180 - ao1
  cao2 = cos(ao2)
  sao2 = sin(ao2)
  
  # Compute the valences {vo1,vo2} of the C-O bonds and {vh1,vh2} of the O-H bonds:
  if eq == 0 :
    vh1 = 1.0  # Valence of lower O-H bond.
    vh2 = 0.0  # Valence of upper O-H bond.
    vo1 = 1.0  # Valence of lower C-O bond.
    vo2 = 2.0  # Valence of upper C-O bond.
  else :
    vh1 = 0.5  # Valence of lower O-H bond.
    vh2 = 0.5  # Valence of upper O-H bond.
    vo1 = 1.5  # Valence of lower C-O bond.
    vo2 = 1.5  # Valence of upper C-O bond.

  # Compute the bond lengths fo the two C-O bonds:
  bo1 = svg.rel_bond_length(vo1)
  bo2 = svg.rel_bond_length(vo2)

  # Compute rise angle {ah1} (radians) and length {bh1} of lower O-H bond:
  if eq == 0 :
    ah1 = AOH*pi/180 - ao1
    bh1 = svg.OHbondlength
  else :
    # Set {bh1} to the ideal length of O-H bond:
    bu = svg.rel_bond_length(1.0)
    bh1 = (0.5*(svg.OHbondlength + bu)/bu)*svg.rel_bond_length(vh1)
    # Compute angle {ah} for mid-position:
    assert math.fabs(ao1-ao2) < 0.0001 # Assume that the oxygens are Y-symmetrical.
    ahmax = 70*pi/180 # Max rise angle of O-H bond.
    ahmin = 45*pi/180 # Min rise angle of O-H bond.
    sah1 = sao1*bo1/bh1
    if sah1 >= 1 :
      ah1 = ahmax
    else :
      ah1 = math.asin(sah1)
      if ah1 > ahmax : ah1 = ahmax
      if ah1 < ahmin : ah1 = ahmin
    # Recompute {bh1} from {bo1} and {ah1}:
    bh1 = sao1*bo1/sin(ah1)
    
  cah1 = cos(ah1)
  sah1 = sin(ah1)

  # Atom centers C,O,O,H:
  pc = [0,0]                                 # Carbon.
  po1 = rn.add(pc,  [+cao1*bo1, -sao1*bo1])  # Lower oxygen.
  po2 = rn.add(pc,  [+cao2*bo2, +sao2*bo2])  # Upper oxygen.
  ph =  rn.add(po1, [+cah1*bh1, +sah1*bh1])  # Hydrogen.

  # Place atoms:
  kc =  fm.add_atom("C", pc,  0,1) # Atom 0.
  ko1 = fm.add_atom("O", po1, 0,1)
  ko2 = fm.add_atom("O", po2, 0,1)
  kh =  fm.add_atom("H", ph,  0,1)
  
  # Insert bonds (with dashes on inner side):
  fm.add_bond(kc,ko1,vo1)
  fm.add_bond(ko2,kc,vo2)
  fm.add_bond(ko1,kh,vh1)
  if vh1 < 0.99 :
    fm.add_bond(kh,ko2,1-vh1)
  
  return fm
  #----------------------------------------------------------------------