from mmtk import *
import string

ff = AmberForceField()
world = InfiniteUniverse(ff)
layer = Group('kaehler_a_layer')
join = Joiner(layer)
join.joinGroups(layer.hookup, layer.hookdown)
join.joinGroups(join.getIntermediateResult().hookup, layer.hookdown)
join.joinGroups(join.getIntermediateResult().hookup, layer.hookdown)
join.joinGroups(join.getIntermediateResult().hookup, layer.hookdown)
result = join.getResult()
try:
    result.setPositions(ff)
except ConformationError, msg:
    print 'ConformationError',msg

world.group = result
minimizer = ConjugateGradientMinimizer(world, log = (0, None, 50, stdout,
						     ("energy",)))
#minimizer(convergence = 1.e-3, steps = 25)
world.removeObject(result)

oxygens = [result.generateNewAtom('O'), result.generateNewAtom('O'),\
	   result.generateNewAtom('O'), result.generateNewAtom('O'),\
	   result.generateNewAtom('O'), result.generateNewAtom('O')]

hpairs = [[result.findAtomByName('kaehler_a_layer.kaehler_a_layer1.H2'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H3')],\
	  [result.findAtomByName('kaehler_a_layer.kaehler_a_layer1.H4'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H5')],\
	  [result.findAtomByName('kaehler_a_layer.kaehler_a_layer1.H6'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H1')],\
	  [result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H6'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer3.H1')],\
	  [result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H4'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer3.H5')],\
	  [result.findAtomByName('kaehler_a_layer.kaehler_a_layer2.H2'),\
	   result.findAtomByName('kaehler_a_layer.kaehler_a_layer3.H3')]]

for hpr in hpairs:
    print hpr[0].fullName(), hpr[1].fullName(), '%.3f' % (hpr[0].position() - hpr[1].position()).length()

result.replaceAtoms(hpairs, oxygens)
for a in oxygens:
    result.setAtomProperty(a, 'amber_atom_type', 'OS')
#for a in result.atoms:
#    print '%3d' % a.index, result.getAtomProperty(a, 'amber_atom_type'), a.fullName(),
#    print hex(id(result.getReference(a))),result.getReference(a),getattr(result, 'amber_atom_type')[result.getReference(a)]

world.group = result
#minimizer = ConjugateGradientMinimizer(world, log = (0, None, 50, stdout,
#						     ("energy",)))
#minimizer(convergence = 1.e-3, steps = 25)

print result.describeBentAngles(ff)

from ConfigIO import OutputFile
fmt = OutputFile('PDB/resultk.pdb', 'pdb.connect')
fmt.write(result)
fmt.close()

for line in result.describeBentAngles(ff):
    print line

if 1:
    import RasMol
    ras_session = RasMol.RasMol(result, 1)
    ras_session.RefreshScreen()
    ras_session.Interact()