#
# This MMTK/Python script is used to generate a pair of "Lego"-block like
# molecules as described in my paper for MNT6.
# writtern by Peter McCluskey (pcm@rahul.net)
#

from mmtk import *
import sys
import Units
import Positioner

ff = GenericForceField('mm3')
world = InfiniteUniverse(ff)
block1 = Group('diamondoid_fragment')
join = Joiner(block1)
width = 5
join.repeatGroup('hookeast', 'hookwest', width)
join.repeatGroup('hooksouth', 'hooknorth', 5, 1, [180*Units.deg]*width)
join.repeatGroup('hookup', 'hookdown', 2)

result = join.getResult()
result.setPositions(ff)

a_hole = result.findAtomByName('diamond4carbon.diamond4carbon7.C4')
surface_point = Positioner.findNearestPointOnSurface(result, a_hole)
surface_vector = surface_point - a_hole.position()
dist_to_surface = surface_vector.length()
surface_dir = surface_vector.normal()
cylinder1 = Cylinder(a_hole.position() - 0.1*Units.Ang*surface_dir,\
		     2.5*Units.Ang, (1*Units.Ang+dist_to_surface)*surface_dir)

atoms_in_cyl1 = result.atomsWithinRegion(cylinder1)
result.removeAtoms(atoms_in_cyl1, None)
c_below = result.findAtomByName('diamond4carbon.diamond4carbon32.C3')
nitrogen = result.generateNewAtom('N')
result.replaceAtoms([c_below], [nitrogen])
result.fillDanglingBonds(ff, 'N', 8, 3)


def filter2n(atom):
    if atom.symbol != 'C':
	return 0
    bonded_to = atom.bondedTo()
    countn = 0
    for a in bonded_to:
	if a.symbol == 'N':
	    countn = countn + 1
    return countn == 2

carbons_bonded_to_2_ns = filter(filter2n , result.atomList())
selected_ns = []
replace_cs = []
for a1 in carbons_bonded_to_2_ns:
    for a2 in a1.bondedTo():
	if a2.symbol == 'C':
	    c = a2
    a1_pos = a1.position()
    for a2 in a1.bondedTo():
	if a2.symbol == 'N':
	    x_product = (a1_pos - c.position()).cross(a2.position() - a1_pos)
	    if surface_dir * x_product > 0:
		selected_ns.append(a2)
    replace_cs.append(result.generateNewAtom('C'))
result.replaceAtoms(selected_ns, replace_cs, ff)
for i in range(len(carbons_bonded_to_2_ns)):        # make double bonds:
    result.attachBondBetween(carbons_bonded_to_2_ns[i], replace_cs[i])

hook3_list = []
for db in result.danglingbonds:
    if nitrogen in db.a1.bondedTo():
	hook3_list.append(db)
join2 = Joiner(result)
cn = Group('cyanide')
for db in hook3_list:
    join2.joinGroups(Hook((result, [db])), cn.hook1, 0)
    result = join2.getIntermediateResult()
result = join2.getResult()
try:
    result.setPositions(ff)
except ConformationError, msg:
    print 'ConformationError',msg

result.hydrogenateDanglingBonds(ff)


result2 = result.clone()	# shouldn't this be Molecule? (no clone yet)
result2.name = 'd42'
face2 = []
for a in result2.atomList():
    if a.symbol == 'N':
	bonded = a.bondedTo()
	if len(bonded) == 1 and bonded[0].symbol == 'C':
	    face2.append(a)

face2 = [face2[0], face2[2], face2[1]]
PositionAdjacent(result,carbons_bonded_to_2_ns, result2,face2, [2*Units.Ang]*3)

from ConfigIO import OutputFile
fmt = OutputFile('PDB/resultl5.1.pdb')
fmt.write(result)
fmt = OutputFile('PDB/resultl5.2.pdb')
fmt.write(result2)
fmt.close()

world.addObject(result)
world.addObject(result2)

fmt = OutputFile('PDB/resultlw.pdb')
fmt.write(world)
fmt.close()

if 0:
    world.view()
else:
    import RasMol
    ras_session = RasMol.RasMol(world, 1)
    ras_session.RefreshScreen()
    ras_session.Interact()