index.html

 avatarjbcooper
lua
11 days ago
47 kB
0
Indexable
Never
<html>
<head><meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>mini map aide - display protein electron density maps on mobile phone browsers - PDB and CCP4 map file viewer for proteins, nucleotides and ligands</title>
<script src='fengari_web.js' type="text/javascript" async></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/three.js/0.155.0/three.min.js"></script>
<script src="http://minimapai.de/surfacenets.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjs/10.1.0/math.min.js"></script>
<style>
  #topbar {background-color:dodgerblue; position:relative; line-height:30px; font-size:22px; text-align:center; color:white;}
  body {margin:auto; font-size:16px; font-family:arial; color:white; background-color:#0C0A0A;}
  p {margin:0px;}
  abbr{border-bottom:none !important; cursor:inherit !important; text-decoration:none !important;}
button:active {color:gray;}
</style></head><body>
<div id="topbar"><p>
<abbr title="Link to help file."><a href="minimaphelp.html" target="info"><button>About Mini Map Aide</button></a></abbr>
<abbr title="The start point in the sequence. This is an internal index for the residues and is not necessarily the same as the residue number in the PDB file."><label for="position" style="margin-left: 10px;">
Go&nbspto: </label><input type="number" min="-999" max="9999" step="1" value="1" id="position" name="position" size="4"></abbr>
<abbr title="The contour level in sigma or rms units. Change sign with +/- button."><label for="contour" style="margin-left: 10px;">
Contour&nbsplevel: <button id="densign" style="margin-left: 10px;">+/-</button></label><input type="number" min="-20" max="20" step="0.1" value="1" id="contour" name="contour" size="4"></abbr>
<abbr title="Select the PDB file."><span style="display:inline-block;"><label for="pdbInput" style="margin-left: 10px;">PDB: </label><input type="file" id="pdbInput" style="width:180px;"></span></abbr>
<abbr title="Select CCP4 map or MTZ file from refinement."><span style="display:inline-block;"><label for="mapInput" style="margin-left: 10px;">Map: </label><input type="file" id="mapInput" style="width:180px;" disabled></span></abbr>
<span style="display:inline-block;"><abbr title="Update display after changing options."><button id="run" style="margin-left: 10px;">Redraw</button></abbr>
<abbr title="Click to view outer part of sidechain or nucleotide base."><button id="viewsidechain" style="margin-left: 10px;">Side chain</button></abbr>
<abbr title="Swap between animation and scrolling."><button id="stop" style="margin-left: 10px;">Scroll</button></abbr>
<abbr title="Move to next residue."><button id="step" style="margin-left: 10px;">Next residue</button></abbr></span></div>
<div id="graphics" style="text-align: center">
<div><span id="info">Select a PDB file</span><span style="position: absolute; right: 0;"><abbr title="Click to move atoms in the residue at the centre of the display. The moveable atoms will then appear green and can be clicked and dragged. The selected atom is enlarged momentarily to indicate which one is moving."><button id="move" style="margin-top:20px; margin-right: 5px;" disabled>Move</button></abbr><abbr title="Tidy up the geometry of the residue after the atoms have been moved. This is a bit approximate and does not use the electron density map at all. It may be repeated if necessary."><button id="tidy" style="margin-right: 5px;" disabled>Tidy</button></abbr><abbr title="Undo the move and tidy operations on this residue."><button id="undo" style="margin-right: 5px;" disabled>Undo</button></abbr><abbr title="Save the changes in a new PDB file. It is *important* to do this before closing the browser tab if any of the changes made during the session are still wanted for future use!"><button id="save" style="margin-right: 5px;" disabled>Save</button></abbr></span></div>
<div id="myCanvas" style="padding: 0; max-width:100%; height:auto;"></div></div>
<script type="application/lua" async>

--  |d|e|f|a|u|l|t|s| |f|o|r| |m|i|n|i| |m|a|p| |a|i|d|e|                       v15sep23

mapcolpos="dodgerblue" mapcolneg="red" bondcol="#FFD633" altbondcol="#0FFF50" maxbondlen=2.0 maxssbondlen=3.0 boxrad=6.5 maprad=4.2

local js = require "js"
local window = js.global
local threejs = window.THREE
local mathjs = window.math
local document = window.document
canvas=document:getElementById("myCanvas")

width=window.innerWidth
height=window.innerHeight
ratio=width/height
renderer=js.new(threejs.WebGLRenderer)
camera=js.new(threejs.PerspectiveCamera,20,ratio,30,50)
--camera.position:set(0,0,10)
camera.position.z=40 
--camp=camera.position print(camp.x,camp.y,camp.z)
scene=js.new(threejs.Scene)
scene:add(camera)
renderer:setSize(width,height)
canvas:appendChild(renderer.domElement)

material=js.new(threejs.LineBasicMaterial)
material.color:set(bondcol)
material.linewidth=2
altmat=js.new(threejs.LineBasicMaterial)
altmat.color:set(altbondcol)
altmat.linewidth=2
matermap=js.new(threejs.LineBasicMaterial)
matermap.linewidth=1

light=js.new(threejs.PointLight,0xFFFFFF)
light.position:set(50, 50, 50);
scene:add(light)

mouseVector = js.new(threejs.Vector2)
projector = js.new(threejs.Raycaster)

mouse_x,mouse_y,theta,offsetX,offsetY=0,0,0,0,0 moving=false spin=true resifound=false viewSideCh=false foundSideCen=false sideAtGrk={"E","D","G","B"}  --ncyc=0 maxdiff=100  --last 2 normally commented out

function mouseDown(event)
MouseDown = true
mouseX = event.offsetX
mouseY = event.offsetY
get_moving_atom()
end

function touchStart(event)
event:preventDefault() -- needed to stop browser reading swipe events
MouseDown = true
mouseX = event.touches[0].pageX - event.touches[0].target.offsetLeft  -- custom event.offsetX,Y needed for mobile
mouseY = event.touches[0].pageY - event.touches[0].target.offsetTop
get_moving_atom()
end

function mouseMove(event)
if MouseDown==true then
deltaX = (event.offsetX - mouseX)
deltaY = (event.offsetY - mouseY)
if moving and intersect[0]~=nil then
	mouse_x=5*(deltaX)/width
	mouse_y=5*(-deltaY)/height
	dist=math.sqrt(mouse_x^2+mouse_y^2)
	if dist>0.1 then update_orig_coords() mouseX=event.offsetX mouseY=event.offsetY end
else
group.rotation.y=(deltaX/300+offsetX) -- swapped deliberately -- don't understand why +offset needed rather than "-" in mobile...
group.rotation.x=(deltaY/300+offsetY)
--	print(group.rotation.y,group.rotation.x)
end
renderer:render(scene, camera)
end
end

function touchMove(event)
if MouseDown==true then
event:preventDefault()  -- needed to stop browser reading swipe events
deltaX = (event.touches[0].pageX - event.touches[0].target.offsetLeft - mouseX)/300  -- custom event.offsetX,Y needed for mobile
deltaY = (event.touches[0].pageY - event.touches[0].target.offsetTop - mouseY)/300
if moving and intersect[0]~=nil then
	mouse_x=500*(deltaX)/width
	mouse_y=500*(-deltaY)/height
	dist=math.sqrt(mouse_x^2+mouse_y^2)
	if dist>0.1 then update_orig_coords() mouseX=event.touches[0].pageX-event.touches[0].target.offsetLeft mouseY=event.touches[0].pageY-event.touches[0].target.offsetTop end
else
group.rotation.y=offsetX-deltaX -- swapped deliberately
group.rotation.x=offsetY+deltaY
end
renderer:render(scene, camera)
end
end

function mouseUp(event)
MouseDown = false
if moving and intersect[0]~=nil then
	update_orig_coords()
else
	offsetX=group.rotation.y -- ditto
	offsetY=group.rotation.x
end
end

function touchEnd(event) mouseUp(event) end


function get_moving_atom()
mouseVector.x = 2 * (mouseX / width) - 1
mouseVector.y = 1 - 2 * ( mouseY / height )
--print(mouseVector.x,mouseVector.y)
if moving then
projector:setFromCamera( mouseVector, camera )
intersect = projector:intersectObject( moveables, true ) --only moveables are clickable
if intersect[0]~=nil then
		--print("intersect object x:",intersect[0].object.position.x)
		movatoms=intersect[0].object
		movatom=movatoms.name --print(movatom)
		--print("movatoms:",window.console:log(movatoms))
		--movatoms.material.color:set("white")
		movatoms.material.size=20
end
end
end


function startmove()
moving=true	
document:getElementById("tidy").disabled=false
document:getElementById("undo").disabled=false
document:getElementById("save").disabled=false
resisafe={} for i,j in pairs(protein[startpos]) do resisafe[i]={} for k,l in pairs(j) do resisafe[i][k]=l end end	
move()
end


function move()
moveables = js.new(threejs.Group)
vecmov=js.new(threejs.Vector3,0,0,0)
geomov=js.new(threejs.BufferGeometry)
for k,l in pairs(protein[startpos]) do
                        matermov=js.new(threejs.PointsMaterial)
                        matermov.size=8
                        matermov.sizeAttenuation=false
                        matermov.color:set("#0FFF50")
			vecmove=window:Array(vecmov:set(l.x-xcent,l.y-ycent,l.z-zcent))
			geomove=geomov:clone()
			geomove:setFromPoints(vecmove)
			meshmove=js.new(threejs.Points,geomove,matermov)
                        meshmove.name=k
			moveables:add(meshmove)
end
group:add(moveables)
renderer:render(scene, camera)
geomov:dispose() geomove:dispose() matermov:dispose() --RAM housekeeping
end


function update_orig_coords()
rotamat=group.matrix  --three.js stores transposed rotation matrix.
rotmat=rotamat.elements
invrot1=window:Array(rotmat[0],rotmat[1],rotmat[2])
invrot2=window:Array(rotmat[4],rotmat[5],rotmat[6])
invrot3=window:Array(rotmat[8],rotmat[9],rotmat[10])
invrot=window:Array(invrot1,invrot2,invrot3)
transprot=js.new(mathjs.matrix,invrot)
mouarr=window:Array(mouse_x,mouse_y,0)
munrot=mathjs:multiply(transprot,mouarr)
munrotx=(mathjs:subset(munrot,mathjs:index(0)))
munroty=(mathjs:subset(munrot,mathjs:index(1)))
munrotz=(mathjs:subset(munrot,mathjs:index(2)))
--print("Unrotated shifts:",munrotx,munroty,munrotz)
protein[startpos][movatom].x=protein[startpos][movatom].x+munrotx
protein[startpos][movatom].y=protein[startpos][movatom].y+munroty
protein[startpos][movatom].z=protein[startpos][movatom].z+munrotz
group:remove(moveables)
move()
end


function save()
format_pdb_normal="%s%5i  %-4s%3s%2s%-4s   %8.3f%8.3f%8.3f%s\n"
format_pdb_metals="%s%5i %-4s %3s%2s%-4s   %8.3f%8.3f%8.3f%s\n"
pdbOutputText=crystline
nAt=1
for i,res in ipairs(protein) do
for j,atom in pairs(res) do
if #j==3 then format_pdb=format_pdb_normal else format_pdb=format_pdb_metals end  --handle left shift of metals
pdbOutputLine=string.format(format_pdb,"ATOM  ",nAt,j,atom.res,atom.chain,atom.num,atom.x,atom.y,atom.z,"  1.00 10.00")
pdbOutputText=pdbOutputText..pdbOutputLine
nAt=nAt+1
end
end
a=document:createElement("a") --create imaginary download button
content=window:Array(pdbOutputText)
contentType=js.new(window.Object)
contentType["type"]="text/plain"
savefile=js.new(window.Blob,content,contentType)
a.href=window.URL:createObjectURL(savefile)
a.download=pdbinput.files[0].name
a:click()
--print("a:",window.console:log(a))
window.URL:revokeObjectURL(a.href)
moving=false
scene:remove(group)
readform()
getBoxAtoms()
getMapBox()
group:add(meshmap)
renderer:render(scene, camera)
document:getElementById("save").disabled=true
end


function undo()
protein[startpos]=resisafe
run()
end


function run()  --redraw
scene:remove(group)
readform()
getBoxAtoms()
if resifound==true then getMapBox() end
document:getElementById("move").disabled=false
document:getElementById("tidy").disabled=true
document:getElementById("undo").disabled=true	
document:getElementById("save").disabled=true	
end


function rerun()  --redraw after tidy
scene:remove(group)
readform()
getBoxAtoms()
if resifound==true then getMapBox() end
document:getElementById("save").disabled=false
end


function step()
scene:remove(group)
document:getElementById("tidy").disabled=true
document:getElementById("undo").disabled=true	
document:getElementById("save").disabled=true
--Commented this block out to avoid problems with insertion codes
--startpos=document:getElementById("position").value
--document:getElementById("position").value=startpos+1
--readform()
startpos=startpos+1
startRes=protein[startpos]
if startRes~=nil then _,firstAtom=next(startRes)
    firstAtomResNum=string.sub(firstAtom.num,1,-2)  firstAtomResNum=tonumber(firstAtomResNum) --remove insertion code so resnum is a number
    document:getElementById("position").value=firstAtomResNum
end
getBoxAtoms()
if resifound==true then getMapBox() end
end


function stop()  --stop/start animation
if spin==true then 
spin=false
renderer:setAnimationLoop() offsetX=group.rotation.y offsetY=group.rotation.x
document:getElementById("stop").innerHTML="Roll"
else 
spin=true
renderer:setAnimationLoop(twirl)
document:getElementById("stop").innerHTML="Scroll"	
end
end


function readfile()
scene:remove(group)
pdbfile=pdbinput.files[0]
--print(pdbfile.name)
reader=js.new(window.FileReader)
--reader.onload=function() print(reader.result) end
reader.onload=readatoms
reader:readAsText(pdbfile)	
end


function readatoms()  --reads all atoms into table protein={res1,res2,..}
nAt=0 oldresnum=0 xtot,ytot,ztot=0,0,0 protein={} residue={} chainold="" chains={} --residue={["CA "]={["res"]=resnam,["num"]=resnum,["x"]=3.1,["y"]=..},["CB "]={},..}
filetext=reader.result
for line in filetext:gmatch("[^\r\n]+") do
        if line:sub(1,6)=="CRYST1" then crystline=line.."\n"
        elseif line:sub(1,6)=="ATOM  " or line:sub(1,6)=="HETATM" then
		atnam=line:sub(14,16)
                altgrp=line:sub(17,17) --print(altgrp)
		attype=atnam:sub(1,1)
                metal=line:sub(13,13) if metal~=" " then atnam=metal..atnam end  --remove left shift of heavier elements
                if attype~="H" and attype~="D" and attype~="G" then  --omit H's and long H-atom names
		--print(line)
		resnam=line:sub(18,20)
                chain=line:sub(21,22)
                resnum=line:sub(23,27)
                if resnum~=oldresnum or chain~=chainold then oldresnum=resnum if nAt>0 then table.insert(protein,residue) residue={} end end
                if chain~=chainold then
                             uniqch=true
                             for i,j in ipairs(chains) do if chain==j then uniqch=false end end
                             if uniqch then table.insert(chains,chain) print("Chain:", chain) end
                             chainold=chain
                end
		x_i=tonumber(line:sub(28,38)) xtot=xtot+x_i
		y_i=tonumber(line:sub(39,46)) ytot=ytot+y_i
		z_i=tonumber(line:sub(47,54)) ztot=ztot+z_i
		atom={["res"]=resnam,["num"]=resnum,["chain"]=chain,["x"]=x_i,["y"]=y_i,["z"]=z_i}
                if altgrp~=" " then atnam=atnam..altgrp end
                residue[atnam]=atom
                nAt=nAt+1
                if nAt==1 then firstAtomResNum=resnum:sub(1,-2)  firstAtomResNum=tonumber(firstAtomResNum) end
		end
	end
end
table.insert(protein,residue) -- applies to last residue only
--print("Number of atoms read: ", n," Last residue number: ",resnum,"\n")	
document:getElementById("mapInput").disabled=false
document:getElementById("info").innerHTML="Select CCP4 map"
chainbox=document:getElementById("position")
chainbox:insertAdjacentHTML("beforebegin","<input value=\""..chains[1].."\" id=\"chain\" name=\"chain\" size=\"2\">")
document:getElementById("position").value=firstAtomResNum
--for i,res in ipairs(protein) do print("Residue:",i) for j,k in pairs(res) do print("Atom:",j) for l,m in pairs(k) do print(l,m) end end end
--for i,res in ipairs(protein) do _,atm=next(res) print(atm.chain,atm.num) end
getBoxAtoms()
end


function getBoxAtoms()  --££
print("Num residues, current index:",#protein,startpos)
if protein[startpos]==nil then document:getElementById("info").innerHTML="Residue not found" resifound=false return
elseif protein[startpos]["CA "]~= nil then centrat=protein[startpos]["CA "] 
else for i,j in pairs(protein[startpos]) do centrat=protein[startpos][i] break end -- if no CA, iterator takes first atom
end
if viewSideCh==true then
for m,n in ipairs(sideAtGrk) do
    for i,j in pairs(protein[startpos]) do
        if i:sub(2,2)==n or i:sub(2,3)=="5 " then  -- go up to 4 steps along a side chain or find atom 5 of base
           centrat=protein[startpos][i]
           foundSideCen=true
           break
        end
    end
if foundSideCen==true then break end
end
end
resifound=true --viewSideCh=false
--print(centrat.x,centrat.y,centrat.z)
xcent,ycent,zcent=centrat.x,centrat.y,centrat.z
xmin,xmax,ymin,ymax,zmin,zmax=centrat.x-boxrad,centrat.x+boxrad,centrat.y-boxrad,centrat.y+boxrad,centrat.z-boxrad,centrat.z+boxrad	

box={}
for i,res in ipairs(protein) do  --##
residue={}
for j,atom in pairs(res) do
	--print(j,atom.x)
	if xmin<=atom.x and atom.x<=xmax and ymin<=atom.y and atom.y<=ymax and zmin<=atom.z and atom.z<=zmax then
		--print("Inside box:",atom.res,atom.num,j)
		residue[j]=atom
	end
end
if next(residue)~=nil then table.insert(box,residue) end		
end		

extrabox={}  -- now only for Cys SG atoms
for i,res in ipairs(box) do  --##  gets the usual bonds
bonds={} cys={} jj=0
for j,atom in pairs(res) do
	jj=jj+1	kk=0 --kk must be reset here (bug found 11aug23)
	if j=="SG " then cys["SG "]=atom table.insert(extrabox,cys) end  -- SG's added to extrabox
	for k,atom2 in pairs(res) do
                other_part=false
                kk=kk+1
		if kk>jj then
                        --print("j=",j,"k=",k,"jj,kk",jj,kk)
                        --if #k==4 and #j==4 then print("kj=",k[4],j[4]) if string.sub(k,4)~=string.sub(j,4) then other_part=true end end
                        if #k==4 and #j==4 then if string.sub(k,4)~=string.sub(j,4) then other_part=true end end
                        d=math.sqrt((atom.x-atom2.x)^2+(atom.y-atom2.y)^2+(atom.z-atom2.z)^2)
			if not other_part and d<maxbondlen then
			bond={{atom.x,atom.y,atom.z},{atom2.x,atom2.y,atom2.z}}
			table.insert(bonds,bond)			
			end
		end
	end
end	
	if res["C  "]~=nil and i<#box then if box[i+1]["N  "]~=nil then 
	if math.sqrt((res["C  "].x-box[i+1]["N  "].x)^2+(res["C  "].y-box[i+1]["N  "].y)^2+(res["C  "].z-box[i+1]["N  "].z)^2) < maxbondlen then
	peptide={{res["C  "].x,res["C  "].y,res["C  "].z},{box[i+1]["N  "].x,box[i+1]["N  "].y,box[i+1]["N  "].z}}
	table.insert(bonds,peptide)
	end end
	end
	if res["O3'"]~=nil and i<#box then if box[i+1]["P  "]~=nil then 
	if math.sqrt((res["O3'"].x-box[i+1]["P  "].x)^2+(res["O3'"].y-box[i+1]["P  "].y)^2+(res["O3'"].z-box[i+1]["P  "].z)^2) < maxbondlen then
	phospho={{res["O3'"].x,res["O3'"].y,res["O3'"].z},{box[i+1]["P  "].x,box[i+1]["P  "].y,box[i+1]["P  "].z}}
	table.insert(bonds,phospho)
	end end
	end 
if #bonds>0 then res["bonds"]=bonds	end			
end  --##

group = js.new(threejs.Group)
drawbonds(box)  -- draws the normal aa bonds

nextra=#extrabox   -- gets SS bonds
if nextra>1 then
    extrabonds={}
    for j,atom in ipairs(extrabox) do
	    for m,n in pairs(atom) do at1=m end  --get atom key
	    for k,atom2 in ipairs(extrabox) do
		    if k>j then
			    for o,p in pairs(atom2) do at2=o end  --ditto
			    d=math.sqrt((extrabox[j][at1].x-extrabox[k][at2].x)^2+(extrabox[j][at1].y-extrabox[k][at2].y)^2+(extrabox[j][at1].z-extrabox[k][at2].z)^2)
			    if d<maxssbondlen then
			    bond={{extrabox[j][at1].x,extrabox[j][at1].y,extrabox[j][at1].z},{extrabox[k][at2].x,extrabox[k][at2].y,extrabox[k][at2].z}}
			    table.insert(extrabonds,bond)
			    end
		    end
	    end
    end
    extraboxbond={{["ssbrs"]={},["bonds"]=extrabonds}}  --dummy residue containing SS bonds
    drawbonds(extraboxbond,true)  -- draws any SSBRs
end
--for i,j in ipairs(extraboxbond) do for k,l in pairs(j) do print(k) if k=="bonds" then for m,n in ipairs(l) do print(n[1][1],n[1][2],n[1][3],n[2][1],n[2][2],n[2][3]) end end end end
drawatoms()
end  --££


function drawbonds(abox,alt)
if alt==true then matter=altmat alt=false else matter=material end
vector1=js.new(threejs.Vector3,0,0,0)
vector2=js.new(threejs.Vector3,0,0,0)
geobond = js.new(threejs.BufferGeometry)	
--for i,res in ipairs(protein) do
for i,res in ipairs(abox) do
	if res["bonds"]~=nil then
	for j,bond in ipairs(res["bonds"]) do 
        vector1:set(bond[1][1]-xcent,bond[1][2]-ycent,bond[1][3]-zcent)
        vector2:set(bond[2][1]-xcent,bond[2][2]-ycent,bond[2][3]-zcent)
		vectors=window:Array(vector1,vector2)
		geometry1=geobond:clone()
		geometry1:setFromPoints(vectors)
		mesh1 = js.new(threejs.Line, geometry1, matter)
		group:add(mesh1)
	end
	end
end
--for i,res in ipairs(box) do print("Residue:",i) for j,k in pairs(res) do print(j) end end
end


function drawatoms() 
vectoratom=js.new(threejs.Vector3,0,0,0)
atgeom=js.new(threejs.BufferGeometry)
for i,j in ipairs(box) do 
		for k,l in pairs(j) do
			attype=k:sub(1,1) 
			--if attype=="N" or attype=="O" or attype=="S" or attype=="P" then
                        if attype~="C" and attype~="b" and attype~="s" then  --"b" and "s" are the dummy atom types for the bonds and SSBR's
                        if attype=="N" then atcol="blue" elseif attype=="O" then atcol="red" elseif attype=="S" then atcol="#0FFF50" else atcol="#0FFF50" end
				--print(attype,atcol,l.x-xcent,l.y-ycent,l.z-zcent)
				materat=js.new(threejs.PointsMaterial)
				materat.size=3
				materat.sizeAttenuation=false
				materat.color:set(atcol)
				vectorat=window:Array(vectoratom:set(l.x-xcent,l.y-ycent,l.z-zcent))
			geomat=atgeom:clone()
			geomat:setFromPoints(vectorat)
			meshat=js.new(threejs.Points,geomat,materat)
			group:add(meshat)
			end
		end
end
--Next lines not needed in map version
scene:add(group)
if spin==true then renderer:setAnimationLoop(twirl)	else group.rotation.x=offsetY group.rotation.y=offsetX renderer:setAnimationLoop() renderer:render(scene, camera) end
document:getElementById("move").disabled=false
end


function twirl()
	theta=0.02
	group.rotation.x=group.rotation.x+theta
	group.rotation.y=group.rotation.y+theta
	renderer:render(scene, camera)
end


function uploadMTZ()
document:getElementById("mtztips").innerHTML="Please wait for uploads to complete!" 
--document:getElementById("mapdownload").disabled=true
diffneeded=document:getElementById("diffmap").checked
formData=js.new(window.FormData)
formData:append("pdbfile",pdbfile)
formData:append("mtzfile",mapfile)
formData:append("diffneeded",diffneeded)
fetchObject=js.new(window.FormData)
fetchObject.method="POST"
fetchObject.body=formData
fetchpromise=window:fetch("/cgi-bin/run_gemmi.py",fetchObject)
fetchpromise["then"](fetchpromise, getMapBack)  --override lua's own if/then command
end


function getMapBack()
document:getElementById("mtztips").innerHTML="Please wait for MAP file to download."
if diffneeded~=true then getmap=nameparts[1]..".map" else getmap=nameparts[1].."_diff.map" end
mapURI="/temp/"..getmap
mapdownload=document:createElement("a")
mapdownload.href=mapURI
mapdownload.download=getmap
mapdownload:click()
document:getElementById("mtztips").innerHTML="Please use map file chooser button to load it."
end


function readmapfile()
mapfile=mapinput.files[0]
name1=pdbfile.name name2=mapfile.name print(name1, name2)
nameparts={} for word in string.gmatch(mapfile.name,"[^.]+") do table.insert(nameparts,word) end --print(nameparts[1],nameparts[2])
ref_filetype=string.lower(nameparts[2])
if ref_filetype=="mtz" and nAt>0 then
	--mtztext="Tick for difference map: <input type=\"checkbox\" id=\"diffmap\">&nbsp &nbsp<button id=\"uploadmtz\">Make map from MTZ file</button>&nbsp &nbsp<button id=\"mapdownload\" disabled>Download map file</button> <p> <small><span id=\"mtztips\"></span></small>"
        mtztext="Tick for difference map: <input type=\"checkbox\" id=\"diffmap\">&nbsp &nbsp<button id=\"uploadmtz\">Make map from MTZ file</button><p><small><span id=\"mtztips\"></span></small>"
        document:getElementById("info").innerHTML=mtztext
	mtzin=document:getElementById("uploadmtz")
	mtzin:addEventListener("click", uploadMTZ)
else	
	mapreader=js.new(window.FileReader)
	mapreader.onload=readmap
	mapreader:readAsArrayBuffer(mapfile)
	document:getElementById("info").innerHTML="Please wait"
end
end


function readmap()
map=mapreader.result
viewer=js.new(window.DataView,map)
print("Grid:")
grid={} for g=1,10 do grid[g]=viewer:getInt32((g-1)*4,true) print(grid[g]) end	--up to 40
na=grid[8] nb=grid[9] nc=grid[10] --grid info num div on axis, min and max grid pts grid[4] is mode i.e. 2 for 32 bit floats
print("Cell:")
cell={} for c=1,6 do cell[c]=viewer:getFloat32(40+(c-1)*4,true) print(cell[c]) end	--up to 64
print("Axis order:")
order={} for o=1,3 do order[o]=viewer:getInt32(64+(o-1)*4,true) print(order[o]) end	--2,1,3 z sections (ncode 1), 3,1,2 y sections (ncode3), 1,2,3 for gemmi, 3,2,1 for EM
print("Min, max, mean density:")
mimameden={} for m=1,3 do mimameden[m]=viewer:getFloat32(76+(m-1)*4,true) print(mimameden[m]) end mean=mimameden[3]
print("Symm:")
sym={} for s=1,3 do sym[s]=viewer:getInt32(88+(s-1)*4,true) print(sym[s]) end denoffset=1024+sym[2] -- sym[1]=sg num sym[2] is length of symop string
rms=viewer:getFloat32(216,true) print("rms:",rms)
nlbl=viewer:getInt32(220,true) --print(nlbl)

lblbufarr=map:slice(224,1023)
lblvwr=js.new(window.TextDecoder,"utf-8")
labels=lblvwr:decode(lblbufarr) --print(labels) 

symbufarr=map:slice(1024,denoffset);
symvwr=js.new(window.TextDecoder,"utf-8")
symops=symvwr:decode(symbufarr) --print(symops) 	
	
--den={} for d=1,6 do den[d]=viewer:getFloat32(denoffset+(d-1)*4,true) print(den[d]) end

if order[1]==2 and order[2]==1 and order[3]==3 then print("Old normal CCP4 map") ncode=1
amin=grid[6] numona=grid[2] bmin=grid[5] numonb=grid[1] cmin=grid[7] numonc=grid[3]
elseif order[1]==3 and order[2]==1 and order[3]==2 then print("Old Mono- or tri-clinic CCP4 map") ncode=3
amin=grid[6] numona=grid[2] bmin=grid[7] numonb=grid[3] cmin=grid[5] numonc=grid[1]
elseif order[1]==1 and order[2]==2 and order[3]==3 then print("Gemmi map") ncode="gemmi"
amin=grid[5] numona=grid[1] bmin=grid[6] numonb=grid[2] cmin=grid[7] numonc=grid[3]
elseif order[1]==3 and order[2]==2 and order[3]==1 then print("CCPEM map") ncode="ccpem"
amin=grid[7] numona=grid[3] bmin=grid[6] numonb=grid[2] cmin=grid[5] numonc=grid[1]
else document:getElementById("info").innerHTML="Bad fast medium slow axis order" return
end
fractorth()
end	

function fractorth()
--xmin,xmax,ymin,ymax,zmin,zmax	-- convert to fractional
a,b,c,alpha_deg,beta_deg,gamma_deg=cell[1],cell[2],cell[3],cell[4],cell[5],cell[6]
alpha=math.rad(alpha_deg) beta=math.rad(beta_deg) gamma=math.rad(gamma_deg)
cos_alpha_star=(math.cos(beta)*math.cos(gamma)-math.cos(alpha))/(math.sin(beta)*math.sin(gamma))
sin_alpha_star=math.sqrt(1-cos_alpha_star^2)
--print("Cos(alpha*), sin(alpha*): ", cos_alpha_star, " ", sin_alpha_star,"<br>")
cos_delta=-math.sin(beta)*cos_alpha_star
cos_epsilon=math.sin(beta)*sin_alpha_star
L1=window:Array(a,b*math.cos(gamma),c*math.cos(beta))
L2=window:Array(0,b*math.sin(gamma),-c*math.sin(beta)*cos_alpha_star)
L3=window:Array(0,0,c*math.sin(beta)*sin_alpha_star)
La=window:Array(L1, L2, L3)
L=js.new(mathjs.matrix,La)	
--print(L) 
Linv=mathjs:inv(L)  --fractionalisation matrix
--print(Linv)
getMapBox()
end
	
function getMapBox()	
document:getElementById("info").innerHTML=centrat.res.." "..centrat.num..centrat.chain
if foundSideCen==false then
   if protein[startpos]["CB "]~=nil then mcenat=protein[startpos]["CB "] elseif protein[startpos]["CA "]~=nil then mcenat=protein[startpos]["CA "] else mcenat=centrat end
else
   mcenat=centrat
   foundSideCen=false
end
print("Map centre:",mcenat.x,mcenat.y,mcenat.z)
xmin=mcenat.x-maprad xmax=mcenat.x+maprad ymin=mcenat.y-maprad ymax=mcenat.y+maprad zmin=mcenat.z-maprad zmax=mcenat.z+maprad
	
clo=window:Array(xmin,ymin,zmin)	
far=window:Array(xmax,ymax,zmax)
fclo=mathjs:multiply(Linv,clo)
ffar=mathjs:multiply(Linv,far)
--print(fclo,ffar)	

oclo=mathjs:multiply(L,fclo)
ofar=mathjs:multiply(L,ffar)
--print(oclo,ofar)

xminf=(mathjs:subset(fclo,mathjs:index(0)))	yminf=(mathjs:subset(fclo,mathjs:index(1))) zminf=(mathjs:subset(fclo,mathjs:index(2)))	
xmaxf=(mathjs:subset(ffar,mathjs:index(0)))	ymaxf=(mathjs:subset(ffar,mathjs:index(1))) zmaxf=(mathjs:subset(ffar,mathjs:index(2)))	
--print(xminf,xmaxf,yminf,ymaxf,zminf,zmaxf)
	
mapamin=math.ceil(xminf*na) intamin=mapamin-amin mapamax=math.ceil(xmaxf*na) intamax=mapamax-amin ndiva=(mapamax-mapamin) --math.abs maybe needed!! for -ve on axis
mapbmin=math.ceil(yminf*nb) intbmin=mapbmin-bmin mapbmax=math.ceil(ymaxf*nb) intbmax=mapbmax-bmin ndivb=(mapbmax-mapbmin)
mapcmin=math.ceil(zminf*nc) intcmin=mapcmin-cmin mapcmax=math.ceil(zmaxf*nc) intcmax=mapcmax-cmin ndivc=(mapcmax-mapcmin)
--print(intamin,intamax,intbmin,intbmax,intcmin,intcmax,"ndivs",ndiva,ndivb,ndivc)

minimap={}	
for  cc=intcmin,intcmax do minimapsection={} for aa=intamin,intamax do minimaprow={} for bb=intbmin,intbmax do
				if ncode==1 then
				posi=denoffset+cc*numona*numonb*4+aa*numonb*4+bb*4
				elseif ncode==3 then
				posi=denoffset+bb*numona*numonc*4+aa*numonc*4+cc*4
				elseif ncode=="gemmi" then
				posi=denoffset+cc*numona*numonb*4+bb*numona*4+aa*4
                                elseif ncode=="ccpem" then
				posi=denoffset+aa*numonc*numonb*4+bb*numonc*4+cc*4
				end
				dens=viewer:getFloat32(posi,true)
				normdens=(dens-mean)/rms
				table.insert(minimaprow,normdens)
end table.insert(minimapsection,minimaprow) end table.insert(minimap,minimapsection) end

--for i,j in ipairs(minimap) do for k,l in ipairs(j) do for m,n in ipairs(l) do print(n)	end end end 
	
ndiv=window:Array(ndiva,ndivb,ndivc)	
start=window:Array(0,0,0)
bounds=window:Array(start,ndiv)
isomesh=window:surfaceNets(ndiv, function(_,x,y,z) return denread(x,y,z) end, bounds)
window.console:log(isomesh)
posis=isomesh.positions	

vertices={}	
	
for i=0,#posis-1 do
	xcontf=(mapamin+posis[i][0])/na  ycontf=(mapbmin+posis[i][1])/nb  zcontf=(mapcmin+posis[i][2])/nc 
	minimapxyzf=window:Array(xcontf,ycontf,zcontf)
	minimapxyzo=mathjs:multiply(L,minimapxyzf)
	minixo=mathjs:subset(minimapxyzo,mathjs:index(0)) miniyo=mathjs:subset(minimapxyzo,mathjs:index(1)) minizo=mathjs:subset(minimapxyzo,mathjs:index(2))
	table.insert(vertices,{minixo-xcent,miniyo-ycent,minizo-zcent})
end

--scene:add(group)

vector1=js.new(threejs.Vector3,0,0,0)
vector2=js.new(threejs.Vector3,0,0,0)
vector3=js.new(threejs.Vector3,0,0,0)
geomap = js.new(threejs.BufferGeometry)

triangs=isomesh.cells		
for i=0,#triangs-1 do
	verts=triangs[i]
	vert1=verts[0]+1 vert2=verts[1]+1 vert3=verts[2]+1  -- +1 to convert to lua indices in vertices table	
	--print("\n",vert1,vert2,vert3)
	vector1:set(vertices[vert1][1],vertices[vert1][2],vertices[vert1][3])
	vector2:set(vertices[vert2][1],vertices[vert2][2],vertices[vert2][3])
	vector3:set(vertices[vert3][1],vertices[vert3][2],vertices[vert3][3])
	vectors=window:Array(vector1,vector2,vector3)
	mapgeo=geomap:clone()
	mapgeo:setFromPoints(vectors)
	meshmap = js.new(threejs.LineLoop, mapgeo, matermap)
	group:add(meshmap)	
end	
scene:add(group)
--if spin==true then renderer:setAnimationLoop(function() twirl() end)	else offsetX,offsetY=0,0 renderer:setAnimationLoop() renderer:render(scene, camera) document:getElementById("move").disabled=false end
if spin==true then renderer:setAnimationLoop(twirl)	else group.rotation.x=offsetY group.rotation.y=offsetX renderer:setAnimationLoop() renderer:render(scene, camera) end
document:getElementById("move").disabled=false

geobond:dispose() geometry1:dispose() atgeom:dispose() geomat:dispose() geomap:dispose() --RAM housekeeping
material:dispose() altmat:dispose() matermap:dispose() matter:dispose() materat:dispose()

end


--Dictionaries
dictable={   --target distances are read from dist table and are grouped to be within about 0.1 Angstroms. *prev/next*=ignored restraints
["N  "]={["CA "]=4,["CB "]=9,["C  "]=9}, --*prev*["C  "],*prev*["CA "],*prev*["O  "]
["CA "]={["N  "]=4,["C  "]=4,["CB "]=4,["O  "]=8,["CG "]=9,["OG "]=8,["OG1"]=8,["CG1"]=9,["CG2"]=9,["SG "]=10}, --*prev*["C  "],*next*["N  "]
["C  "]={["CA "]=4,["CB "]=9,["N  "]=9,["O  "]=1}, --*next*["N  "],*next*["CA "],*next pro*["CD "]
["O  "]={["C  "]=1,["CA "]=8}, --*next*["N  "]
["CB "]={["CA "]=4,["N  "]=9,["C  "]=9,["CG "]=4,["OG "]=3,["OG1"]=3,["CG1"]=4,["CG2"]=4,["SG "]=5,["CD "]=9,["CD1"]=9,["CD2"]=9,["OD1"]=8,["OD2"]=8,["ND2"]=8,["ND1"]=9,["SD "]=10},
["CG "]={["CA "]=9,["CB "]=4,["CD1"]=4,["CD2"]=4,["CD "]=4,["SD "]=5,["OD1"]=1,["OD2"]=1,["ND2"]=2,["ND1"]=2,["OE1"]=8,["OE2"]=8,["NE1"]=6,["NE2"]=8,["NE "]=9,["CE1"]=8,["CE2"]=8,["CE3"]=9,["CE "]=10}, --Met/Lys CG/CE sorted in code
["OG "]={["CB "]=3,["CA "]=8},
["SG "]={["CB "]=5,["CA "]=10},
["OG1"]={["CB "]=3,["CA "]=8,["CG2"]=8},
["CG1"]={["CB "]=4,["CD1"]=4,["CA "]=9,["CG2"]=9},
["CG2"]={["CB "]=4,["CA "]=9,["OG1"]=8,["CG1"]=9},
["CD "]={["CG "]=4,["CB "]=9,["OE1"]=1,["OE2"]=1,["NE2"]=2,["CE "]=4,["NE "]=4,["CZ "]=9,["NZ "]=9}, --*pro prev*["C  "]
["CD1"]={["CG "]=4,["CG1"]=4,["CB "]=9,["CZ "]=8,["CD2"]=8,["CE1"]=3,["NE1"]=3},   --cd1/2 aliphatic/aromatic simplified
["CD2"]={["CG "]=4,["CE2"]=3,["CB "]=9,["CZ "]=8,["CE3"]=3,["NE1"]=6,["CZ2"]=9,["CZ3"]=8,["NE2"]=3,["ND1"]=6,["CD1"]=8},
["OD1"]={["CG "]=1,["CB "]=8,["OD2"]=6,["ND2"]=6},
["OD2"]={["CG "]=1,["CB "]=8,["OD1"]=6},
["ND1"]={["CG "]=3,["CB "]=9,["CE1"]=2,["CD2"]=6,["NE2"]=6},
["ND2"]={["CG "]=2,["CB "]=8,["OD1"]=6},
["SD "]={["CG "]=5,["CE "]=5,["CB "]=10},
["CE "]={["SD "]=5,["NZ "]=4,["CD "]=4,["CG "]=10},
["CE1"]={["CD1"]=3,["CG "]=8,["CZ "]=3,["OH "]=8,["ND1"]=2,["NE2"]=2,["CE2"]=8},  --ce1/2 same as cd1/2
["CE2"]={["CD2"]=3,["CG "]=8,["CZ "]=3,["OH "]=8,["NE1"]=3,["CZ2"]=3,["CH2"]=8,["CE3"]=8,["CE1"]=8},
["CE3"]={["CD2"]=3,["CG "]=9,["CZ3"]=3,["CH2"]=8,["CE2"]=8},
["OE1"]={["CD "]=1,["CG "]=8,["OE2"]=6,["NE2"]=6},
["OE2"]={["CD "]=1,["CG "]=8,["OE1"]=6},
["NE "]={["CD "]=4,["CG "]=9,["CZ "]=2,["NH1"]=7,["NH2"]=7},
["NE1"]={["CD1"]=3,["CG "]=6,["CE2"]=3,["CD2"]=6,["CZ2"]=9},
["NE2"]={["CD "]=2,["CG "]=8,["OE1"]=6,["CD2"]=2,["CE1"]=2,["ND1"]=6},
["CZ "]={["CE1"]=3,["CE2"]=3,["CD1"]=8,["CD2"]=8,["OH "]=3,["NE "]=2,["CD "]=9,["NH1"]=2,["NH2"]=2},
["NZ "]={["CE "]=4,["CD "]=9},
["CZ2"]={["CE2"]=3,["CD2"]=9,["NE1"]=9,["CH2"]=3,["CZ3"]=8},
["CZ3"]={["CE3"]=3,["CH2"]=3,["CD2"]=8,["CZ2"]=8},
["CH2"]={["CZ2"]=3,["CZ3"]=3,["CE2"]=8,["CE3"]=8},
["NH1"]={["CZ "]=2,["NH2"]=7,["NE "]=7},
["NH2"]={["CZ "]=2,["NH1"]=7,["NE "]=7},
["OH "]={["CZ "]=3,["CE1"]=8,["CE2"]=8}
}

protable={   --target distances for proline
["N  "]={["CA "]=4,["CB "]=8,["C  "]=9,["CG "]=8,["CD "]=4}, --*prev*["C  "],*prev*["CA "],*prev*["O  "]
["CA "]={["N  "]=4,["C  "]=4,["CB "]=4,["O  "]=8,["CG "]=8,["CD "]=9}, --*prev*["C  "],*next*["N  "]
["C  "]={["CA "]=4,["CB "]=9,["N  "]=9,["O  "]=1}, --*next*["N  "],*next*["CA "],*next pro*["CD "]
["O  "]={["C  "]=1,["CA "]=8}, --*next*["N  "]
["CB "]={["CA "]=4,["N  "]=8,["C  "]=9,["CG "]=4,["CD "]=9},
["CG "]={["CA "]=8,["CB "]=4,["CD "]=4,["N  "]=8},
["CD "]={["CG "]=4,["CB "]=9,["N  "]=4,["CA "]=9} --*pro prev*["C  "]
}

tetrigtable={  --tetragonal or trigonal carbon centres, neighbours and distance from centroid of neighbours. CA and Thr/Ile CB chirality not checked.
["ASP"]={["CG "]={"CB ","OD1","OD2",0.1}},["ASN"]={["CG "]={"CB ","OD1","ND2",0.1}},
["GLU"]={["CD "]={"CG ","OE1","OE2",0.1}},["GLN"]={["CD "]={"CG ","OE1","NE2",0.1}},
["PHE"]={["CG "]={"CB ","CD1","CD2",0.0}},
["TYR"]={["CG "]={"CB ","CD1","CD2",0.0},["CZ "]={"CE1","CE2","OH ",0.0}},
["TRP"]={["CG "]={"CB ","CD1","CD2",0.07},["CD2"]={"CG ","CE2","CE3",0.1},["CE2"]={"CD2","NE1","CZ2",0.1}},
["LEU"]={["CG "]={"CB ","CD1","CD2",0.5}},["VAL"]={["CB "]={"CA ","CG1","CG2",0.5}},
["ILE"]={["CB "]={"CA ","CG1","CG2",0.5}},["THR"]={["CB "]={"CA ","OG1","CG2",0.5}},
["ARG"]={["CZ "]={"NE ","NH1","NH2",0.0}},["HIS"]={["CG "]={"CB ","ND1","CD2",0.04}}}

planetable={  --restrain some interatomic vectors to be parallel e.g. CD1-CE1 bond is parallel to CB-CG in Phe/Tyr but shorter by factor of 0.91.
["PHE"]={["CE1"]={"CB ","CG ","CD1",0.91},["CE2"]={"CB ","CG ","CD2",0.91},["CZ "]={"CG ","CD1","CE2",1.0},["CZ "]={"CG ","CD2","CE1",1.0}},
["TYR"]={["CE1"]={"CB ","CG ","CD1",0.91},["CE2"]={"CB ","CG ","CD2",0.91},["CZ "]={"CG ","CD1","CE2",1.0},["CZ "]={"CG ","CD2","CE1",1.0},["OH "]={"CB ","CG ","CZ ",0.93}},
["TRP"]={["CH2"]={"CD2","CE2","CZ3",1.0},["CH2"]={"CD2","CE3","CZ2",1.0},["CZ3"]={"CE2","CZ2","CE3",1.0},["CZ2"]={"CG ","NE1","CE3",1.26},["CH2"]={"CG ","NE1","CZ3",0.63}},
["ARG"]={["NH2"]={"CD ","NE ","CZ ",1.0}}}


function addTetrig()  --calculate distance of central atom from centroid of neighbours and the diff from ideal
xcentet=(residue[nayb1].x+residue[nayb2].x+residue[nayb3].x)/3
ycentet=(residue[nayb1].y+residue[nayb2].y+residue[nayb3].y)/3
zcentet=(residue[nayb1].z+residue[nayb2].z+residue[nayb3].z)/3
dx=residue[centat].x-xcentet dy=residue[centat].y-ycentet dz=residue[centat].z-zcentet
dcalc=math.sqrt(dx^2+dy^2+dz^2)
diff=dtarg-dcalc
table.insert(atDeltas,{["diff"]=diff,["at1"]="tetrig",["at2"]=centat,["delx"]=diff*dx/dcalc,["dely"]=diff*dy/dcalc,["delz"]=diff*dz/dcalc})
end

function addPlane()  --calculate distance of atom from ideal position that would make the vectors parallel
planx=residue[tail2].x+(residue[tip1].x-residue[tail1].x)*dscale
plany=residue[tail2].y+(residue[tip1].y-residue[tail1].y)*dscale
planz=residue[tail2].z+(residue[tip1].z-residue[tail1].z)*dscale
dx=residue[planat].x-planx dy=residue[planat].y-plany dz=residue[planat].z-planz
dcalc=math.sqrt(dx^2+dy^2+dz^2)
diff=dtarg-dcalc
table.insert(atDeltas,{["diff"]=diff,["at1"]="plane",["at2"]=planat,["delx"]=diff*dx/dcalc,["dely"]=diff*dy/dcalc,["delz"]=diff*dz/dcalc})
end

function addTrpPlane()  --calculate distance of CD1 from its ideal position on the line bisecting CH2-CZ3 and CD2-CE2 bonds. Helps to make the Trp flat.
dum1x=(residue["CZ3"].x+residue["CH2"].x)/2 dum1y=(residue["CZ3"].y+residue["CH2"].y)/2 dum1z=(residue["CZ3"].z+residue["CH2"].z)/2
dum2x=(residue["CD2"].x+residue["CE2"].x)/2 dum2y=(residue["CD2"].y+residue["CE2"].y)/2 dum2z=(residue["CD2"].z+residue["CE2"].z)/2
planx=dum1x+(dum2x-dum1x)*1.89 plany=dum1y+(dum2y-dum1y)*1.89 planz=dum1z+(dum2z-dum1z)*1.89
dx=residue["CD1"].x-planx dy=residue["CD1"].y-plany dz=residue["CD1"].z-planz
dcalc=math.sqrt(dx^2+dy^2+dz^2)
diff=dtarg-dcalc
table.insert(atDeltas,{["diff"]=diff,["at1"]="plane",["at2"]="CD1",["delx"]=diff*dx/dcalc,["dely"]=diff*dy/dcalc,["delz"]=diff*dz/dcalc})
end

function checkChir()  --calculate the chiral volume
print(chirat,nayb1,nayb2,nayb3)
chrat=residue[chirat] nbr1=residue[nayb1]  nbr2=residue[nayb2]  nbr3=residue[nayb3]
detRow1=window:Array(nbr1.x-chrat.x, nbr1.y-chrat.y, nbr1.z-chrat.z)
detRow2=window:Array(nbr2.x-chrat.x, nbr2.y-chrat.y, nbr2.z-chrat.z)
detRow3=window:Array(nbr3.x-chrat.x, nbr3.y-chrat.y, nbr3.z-chrat.z)
chirDet=window:Array(detRow1,detRow2,detRow3)
chivdet=js.new(mathjs.matrix,chirDet)
chiv=mathjs:det(chivdet)
print("Chiv:",chiv)
if chiv<0 then correctChir() end
end

function correctChir()
xcentet=(nbr1.x+nbr2.x+nbr3.x)/3  ycentet=(nbr1.y+nbr2.y+nbr3.y)/3  zcentet=(nbr1.z+nbr2.z+nbr3.z)/3
mov1=chrat  mov2=nbr3
dx=mov1.x-xcentet dy=mov1.y-ycentet dz=mov1.z-zcentet  dcalc=math.sqrt(dx^2+dy^2+dz^2)  diff=-0.5 --diff=dtarg-dcalc
delx=diff*dx/dcalc  dely=diff*dy/dcalc  delz=diff*dz/dcalc
mov1.x=mov1.x+delx  mov1.y=mov1.y+dely  mov1.z=mov1.z+delz
mov2.x=mov2.x-delx  mov2.y=mov2.y-dely  mov2.z=mov2.z-delz
end


function mmm()    --massively minimal minimiser (mmm) mainly for side chains (pat not pending)  

if ncyc>0 and ncyc%20==0 and resnam~="GLY" then chirat="CA " nayb1="N  " nayb2="C  " nayb3="CB " checkChir() end  --Check the main chain chiral volume at intervals
if resnam=="THR" or resnam=="ILE" then  --Check the side chain chiral volume at bigger intervals
  if ncyc>0 and ncyc%33==0 then
    chirAtObjs=tetrigtable[resnam]
    for k,l in pairs(chirAtObjs) do  --$$
        --print(k,l)
        chirat=k nayb1=l[1] nayb2=l[2] nayb3=l[3] checkChir()
    end  --$$
  end
end

dist={1.2,1.3,1.4,1.5,1.8,2.2,2.3,2.4,2.5,2.8}
--     1   2   3   4   5   6   7   8   9  10   simplified distances
atDeltas={}
--for k,l in ipairs(lookup1) do print(l) end
kk=0
for k,l in pairs(residue) do
kk=kk+1  mm=0
for m,n in pairs(residue) do
mm=mm+1

if mm>kk then   --calculate the deviations from target distances (diff)
      --print(k,l.x,l.y,l.z)
      if dictable[k]~=nil and dictable[m]~=nil then
          if resnam~="PRO" then dtarg=dist[dictable[k][m]] else dtarg=dist[protable[k][m]] end
          if resnam=="LYS" then
					    if k=="CG " and m=="CE " then dtarg=2.5
					    elseif k=="CE " and m=="CG " then dtarg=2.5
					    end
          end
          if dtarg~=nil then  --print(k,m)
                 dx=n.x-l.x
                 dy=n.y-l.y
                 dz=n.z-l.z
                 dcalc=math.sqrt(dx^2+dy^2+dz^2)
                 diff=dtarg-dcalc
                 table.insert(atDeltas,{["diff"]=diff,["at1"]=k,["at2"]=m})
          end
      end
end

end
end

--Calculate the tetragonal and trigonal target distances
if resnam~="GLY" then centat="CA " nayb1="N  " nayb2="CB " nayb3="C  " dtarg=0.5 addTetrig() end
if tetrigtable[resnam]~=nil then  --@@
    centAtObjs=tetrigtable[resnam]
    for k,l in pairs(centAtObjs) do  --$$
        centat=k nayb1=l[1] nayb2=l[2] nayb3=l[3] dtarg=l[4] addTetrig()
    end  --$$
end  --@@

if ncyc>0 and maxdiff<0.3 and ncyc%3==0 then  --trickle in the planarity targets
if planetable[resnam]~=nil then  --@@
    planeAtObjs=planetable[resnam]
    for k,l in pairs(planeAtObjs) do  --$$
        planat=k tail1=l[1] tip1=l[2] tail2=l[3] dscale=l[4] addPlane()
    end  --$$
end  --@@
if resnam=="TRP" then addTrpPlane() end --extra planarity target for CD1 of Trp
end

ncyc=ncyc+1
--print("Cycle number:",ncyc)
--for i,j in ipairs(atDeltas) do for k,l in pairs(j) do print(l) end end
maxdiff=-100
for i,j in ipairs(atDeltas) do absdiff=math.abs(j.diff)  --find the worst distance violation
    if absdiff>maxdiff then
        maxdiff=absdiff
        worstPair={j.at1,j.at2}
        if j.at1=="tetrig" or j.at1=="plane" then delx=j.delx dely=j.dely delz=j.delz
        end
    end
end

if ncyc>=100 and maxdiff<0.1 then sorted=true end
--for i,j in ipairs(worstPair) do print(i,j) end
at1=worstPair[1]
at2=worstPair[2]
--print("The worst pair (atom name) with delta:",at1,at2,maxdiff)
if at1=="tetrig" or at1=="plane" then --##
mov=residue[at2]
--print(at1)
else
--print("d")
sporad1=math.random(1,2)  --decide at random which atom to move for atom-atom distance violations.
if sporad1==1 then fixed=at1 moving=at2 else fixed=at2 moving=at1 end
if moving=="N  " then moving=fixed fixed="N  " end  --main chain N and C fixed
if moving=="C  " then moving=fixed fixed="C  " end
--print("final decision fixed:",fixed,"moving:",moving)
fix=residue[fixed]
mov=residue[moving]
dx=mov.x-fix.x
dy=mov.y-fix.y
dz=mov.z-fix.z
dcalc=math.sqrt(dx^2+dy^2+dz^2)
if resnam~="PRO" then dtarg=dist[dictable[fixed][moving]] else dtarg=dist[protable[fixed][moving]] end
if resnam=="LYS" then 
	if fixed=="CG " and moving=="CE " then dtarg=2.5
	elseif fixed=="CE " and moving=="CG " then dtarg=2.5
	end
end		
diff=dtarg-dcalc
delx=diff*dx/dcalc dely=diff*dy/dcalc delz=diff*dz/dcalc
end  --##

mov.x=mov.x+delx/2  mov.y=mov.y+dely/2  mov.z=mov.z+delz/2  -- Dampen the shifts
--mov.x=mov.x+delx  mov.y=mov.y+dely  mov.z=mov.z+delz

end  --mmm


function tidy()
residue=protein[startpos]
resnam=string.sub(residue["CA "].res,1,3)
ncyc=0   --normally needed
sorted=false
--mmm() --normally commented out
--repeat mmm() until sorted or ncyc>400
repeat mmm() until sorted or ncyc>200
print("Residue:",resnam,startpos,"max delta:",maxdiff,"after", ncyc, "cycles.")
rerun() --needed normally
end


function denread(xiso,yiso,ziso)
--print("denread:",xiso,yiso,ziso)	
return minimap[ziso+1][xiso+1][yiso+1]-contour_level		
end	


function getStartPos()  --if startpos < 1st resnum go to 1st res in chain, if startpos > last resnum go to last res in chain
for i,res in ipairs(protein) do
    _,firstAtom=next(res) --print(firstAtom.chain,firstAtom.num)
    firstAtomResNum=string.sub(firstAtom.num,1,-2)  firstAtomResNum=tonumber(firstAtomResNum) --remove insertion code so resnum is a number
    if chain==firstAtom.chain then
       lastInChain=i
       if startpos<=firstAtomResNum then wantedpos=i break else wantedpos=lastInChain end
    end
end
startpos=wantedpos  print("Actual res num, index:",firstAtomResNum,wantedpos)
document:getElementById("position").value=firstAtomResNum  --send resnum to the form ignoring insertion code
end


function readform()
document:getElementById("info").innerHTML="Please wait"
position=document:getElementById("position")
startpos=tonumber(position.value)
chain=document:getElementById("chain")
if chain~=js.null then chain=chain.value chain=string.upper(chain) chain=string.format("%2s",chain) getStartPos() end
contour=document:getElementById("contour")
contour_level=tonumber(contour.value) if contour_level<0 then mapcol=mapcolneg else mapcol=mapcolpos end 
matermap.color:set(mapcol)
moving=false
end


function swapdensign()
contour_level=tonumber(document:getElementById("contour").value)
contour_level=-contour_level
document:getElementById("contour").value=contour_level
--print(document:getElementById("contour").value)
end


sidechbutton=document:getElementById("viewsidechain")
sidechbutton:addEventListener("click", function() if viewSideCh==false then sidechbutton.innerHTML="Main chain" viewSideCh=true run() else viewSideCh=false sidechbutton.innerHTML="Side chain" run() end end)

runbutton=document:getElementById("run")
runbutton:addEventListener("click", run)

stepbutton=document:getElementById("step")
stepbutton:addEventListener("click", step)

stopbutton=document:getElementById("stop")
stopbutton:addEventListener("click", stop)

movebutton=document:getElementById("move")
movebutton:addEventListener("click", startmove)

savebutton=document:getElementById("tidy")
savebutton:addEventListener("click", tidy)

undobutton=document:getElementById("undo")
undobutton:addEventListener("click", undo)

savebutton=document:getElementById("save")
savebutton:addEventListener("click", save)

pdbinput=document:getElementById("pdbInput")
pdbinput:addEventListener("change", function() readform() readfile() end)

mapinput=document:getElementById("mapInput")
mapinput:addEventListener("change", function() readform() readmapfile() end)

posneg=document:getElementById("densign")
posneg:addEventListener("click", swapdensign)

canvas:addEventListener("mousedown", function(_,event) mouseDown(event) end)
canvas:addEventListener("mousemove", function(_,event) mouseMove(event) end)
canvas:addEventListener("mouseup",   function(_,event) mouseUp(event)   end)

canvas:addEventListener("touchstart", function(_,event) touchStart(event) end) -- add ",false" after "end" to stop browser reading swipe events?
canvas:addEventListener("touchmove", function(_,event) touchMove(event) end) -- ditto
canvas:addEventListener("touchend",  function(_,event) touchEnd(event)   end) -- ditto

</script>
  </body>
</html>