! 'Hingefind', partitioning algorithm with rotational fitting
! This software is copyrighted, (c) 1995, by Willy Wriggers under the
! terms of the legal statement in the distribution.
! Version 6-22-95.

set echo=off end
set messages=off end

evaluate ($ndomains = 5)	! number of domains to be assigned
evaluate ($maxccounter = 20)	! number of max iterations in partitioning
evaluate ($assign = "slo")	! if "man", upto 9 domains assigned manually
				! if "fas", fast partitioning 
				! if "slo", slow part. with connected domains 
evaluate ($nndist = 2)		! max allowed distance of connected residues
				! for "slo" option, should not be changed

! load topology & parameter files, structure:

@partition.str

! assign manually (gets ignored if $assign # "man"):

vector identify (store1) ((segid AP0)and(resname ALA)) 
vector identify (store2) ((segid AP0)and(resname VAL))
vector identify (store3) ((segid AP0)and(resname THR))
vector identify (store4) (not(ALL))
vector identify (store5) (not(ALL))
vector identify (store6) (not(ALL))
vector identify (store7) (not(ALL))
vector identify (store8) (not(ALL))
vector identify (store9) (not(ALL))

! load resolution of partitioning

@casefile		! percentage of tolerance, rel. initial rmsd

! load coordinates

@casefile1
@casefile2
evaluate ($case1COO = "/yourpath/" + $case1 + 
        "/" + $case1 + ".COO")
coordinates
	@$case1COO	
evaluate ($case2COO = "/yourpath/" + $case2 + 
        "/" + $case2 + ".COO")
coordinates
	disposition = comparison
	@$case2COO	
evaluate ($fname = encode ($tolperc))
evaluate ($oname = "/yourpath/Partition/" + "hinge" + $fname + "." + $case1 + "." + $case2 + ".pdb")
evaluate ($uname = "/yourpath/Partition/" + "hinge" + $fname + "." + $case1 + "." + $case2 + ".psf")
evaluate ($dname = "/yourpath/Partition/" + "hinge" + $fname + "." + $case1 + "." + $case2 + ".dat")
set display $dname end		! output file

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! partitioning algorithm starts here
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

! display some output

display -----------------------------------------------------
display         H   I   N   G   E   F   I   N   D
display ---------------- XPLOR 3.1, $SYSTEM -------------------
display 
display Filename: 
display $dname 
display Created $DATE by user $NAME at $TIME  
display Compares the two cases $case1 and $case2  
display Corresponding filenames of cases:
display $case1COO
display $case2COO 
display Number of domains = $ndomains 

if ($assign = "man") then

  if ($ndomains GE 1) then vector do (segid = AP1) (store1) end if
  if ($ndomains GE 2) then vector do (segid = AP2) (store2) end if
  if ($ndomains GE 3) then vector do (segid = AP3) (store3) end if
  if ($ndomains GE 4) then vector do (segid = AP4) (store4) end if
  if ($ndomains GE 5) then vector do (segid = AP5) (store5) end if
  if ($ndomains GE 6) then vector do (segid = AP6) (store6) end if
  if ($ndomains GE 7) then vector do (segid = AP7) (store7) end if
  if ($ndomains GE 8) then vector do (segid = AP8) (store8) end if
  if ($ndomains = 9) then vector do (segid = AP9) (store9) end if
  coordinates
	fit	
	selection = (all)
	mass = true
	lsq = true
  end
  coordinates
	rms
	mass = true
  end
  evaluate ($thisrmsd = $RESULT)
  display Initial rmsd between cases = $thisrmsd Angstrom.
  display Manual assignment of domains.
  coordinates
	fit	
	selection = (segid AP1)
	mass = true
	lsq = true
  end


elseif ($assign = "slo")  then

  vector identify (store1) (segid AP0)
  coordinates
	fit	
	selection = (store1)
	mass = true
	lsq = true
  end
  coordinates
	rms
	mass = true
	selection = (store1)
  end
  evaluate ($thisrmsd = $RESULT)
  evaluate ($tolerance = $tolperc * $thisrmsd / 100)
  display Initial rmsd between cases = $thisrmsd Angstrom
  display Automatic partitioning, resolution = $tolerance A ( $tolperc % )
  display Slow algorithm, connectivity of domains maintained.
  evaluate ($domain = 1)
  vector identify (store2) (store1)
  vector identify (store5) (not(all))
  vector do (QCOMP = 999) (all)	! old rmsd of converged set large initially

  while ($domain LE $ndomains) loop newdomain
	
	vector identify (store2) ((store2)and(not(store5)))
	evaluate ($num = encode($domain))
	evaluate ($sena = "AP" + $num)
	coordinates
		fit	
		selection = (store2)
		mass = true
		lsq = true
	end
	coordinates
		rms
		mass = true
		selection = (store1)
	end
	evaluate ($status = "converging")
	evaluate ($ccounter = 1)
	while ($status # "converged") loop converge
		vector identify (store9) (tag and (store2))
		for $atom_id in id (store9) loop avres
			vector show norm (RMSD) (byresidue (id $atom_id))
			vector do (BCOMP = $RESULT) (byresidue (id $atom_id))
                end loop avres
  		vector identify (store3) 
				((store2)and(attribute BCOMP < $tolerance))
		vector identify (store5) (not(all))
		vector show sum (mass) (store3)
		evaluate ($mass3 = $RESULT)
		evaluate ($mass5 = 0)
		while ($mass3 > $mass5) loop comparemass
			vector show element (resid) (tag and (store3))
			vector identify (store4) ((resid $RESULT)and(store3))
			vector show sum (mass) (store4)
			evaluate ($mass4 = $RESULT)
			evaluate ($massold = 0)
			while ($mass4 # $massold) loop grow
				evaluate ($massold = $mass4)
				vector identify (store4)
			 		((byresidue((store4) around
					$nndist))and(store3))
				vector show sum (mass) (store4)
				evaluate ($mass4 = $RESULT) 
			end loop grow
			if ($mass4 > $mass5) then 
				vector identify (store5) (store4)
				evaluate ($mass5 = $mass4)	
			end if
			vector identify (store3) ((store3)and(not(store4)))
			evaluate ($mass3 = $mass3 - $mass4)
		end loop comparemass
		coordinates
			fit	
			selection = (store5)
			mass = true
			lsq = true
		end
			coordinates
			rms
			mass = true
			selection = (store1)
		end
		evaluate ($conv1 = 
			$THETA1*$THETA1+$THETA2*$THETA2+$THETA3*$THETA3)
 		evaluate ($conv2 = ($X*$X+$Y*$Y+$Z*$Z)*10000 + $conv1)		
 		if ($conv2 < 100) then evaluate ($status = "converged")	
		end if
		if ($ccounter = $maxccounter) then 
			evaluate ($status = "converged")
			display
			display WARNING: segid $sena not converged!	
			display
		end if
		evaluate ($ccounter = $ccounter + 1)

	end loop converge
{}	
	! label all residues in current domain which meet tolerance criteria
	! and have smaller rmsd than previous domains
	! (optional, increases size of domains AP2, AP3,... !)
		
	vector identify (store8) (tag and (segid AP*))
	for $atom_id in id (store8) loop avres
		vector show norm (RMSD) (byresidue (id $atom_id))
		vector do (BCOMP = $RESULT) (byresidue (id $atom_id))
	end loop avres
	vector do (HARM = BCOMP - QCOMP) (all)
	evaluate ($massold = 0)
	while ($mass5 # $massold) loop grow
		evaluate ($massold = $mass5)
		vector identify (store5)
			((byresidue((store5) around $nndist))
			and((attr BCOMP < $tolerance)and(attr HARM < 0)))
		vector show sum (mass) (store5)
		evaluate ($mass5 = $RESULT) 
	end loop grow
	vector do (QCOMP = BCOMP) (store5)
{}        
        vector do (segid = $sena) (store5)
        evaluate ($domain = $domain + 1)

  end loop newdomain

else

  vector identify (store1) (segid AP0)
  coordinates
	fit	
	selection = (store1)
	mass = true
	lsq = true
  end
  coordinates
	rms
	mass = true
	selection = (store1)
  end
  evaluate ($thisrmsd = $RESULT)
  evaluate ($tolerance = $tolperc * $thisrmsd / 100)
  display Initial rmsd between cases = $thisrmsd Angstrom
  display Automatic partitioning, resolution = $tolerance A ( $tolperc % )
  display Fast algorithm, connectivity of domains not maintained.
  evaluate ($domain = 1)
  vector identify (store2) (store1)
  vector identify (store5) (not(all))
  vector do (QCOMP = 999) (all)	! old rmsd of converged set large initially

  while ($domain LE $ndomains) loop newdomain
	
	vector identify (store2) ((store2)and(not(store5)))
	evaluate ($num = encode($domain))
	evaluate ($sena = "AP" + $num)
	coordinates
		fit	
		selection = (store2)
		mass = true
		lsq = true
	end
	coordinates
		rms
		mass = true
		selection = (store1)
	end
	evaluate ($status = "converging")
	evaluate ($ccounter = 1)
	while ($status # "converged") loop converge
		vector identify (store9) (tag and (store2))
		for $atom_id in id (store9) loop avres
			vector show norm (RMSD) (byresidue (id $atom_id))
			vector do (BCOMP = $RESULT) (byresidue (id $atom_id))
                end loop avres
  		vector identify (store5) 
				((store2)and(attribute BCOMP < $tolerance))
 		coordinates
			fit	
			selection = (store5)
			mass = true
			lsq = true
		end
			coordinates
			rms
			mass = true
			selection = (store1)
		end

		evaluate ($conv1 = 
			$THETA1*$THETA1+$THETA2*$THETA2+$THETA3*$THETA3)
 		evaluate ($conv2 = ($X*$X+$Y*$Y+$Z*$Z)*10000 + $conv1)		
 		if ($conv2 < 100) then evaluate ($status = "converged")	
		end if
		if ($ccounter = $maxccounter) then 
			evaluate ($status = "converged")
			display
			display WARNING: segid $sena not converged!	
			display
		end if
		evaluate ($ccounter = $ccounter + 1)

	end loop converge
{}
	! label all residues in current domain which meet tolerance criteria
	! and have smaller rmsd than previous domains
	! (optional, increases size of domains AP2, AP3,... !)
		
	vector identify (store8) (tag and (segid AP*))
	for $atom_id in id (store8) loop avres
		vector show norm (RMSD) (byresidue (id $atom_id))
		vector do (BCOMP = $RESULT) (byresidue (id $atom_id))
	end loop avres
	vector do (HARM = BCOMP - QCOMP) (all)
  	vector identify (store5) 
		((attr BCOMP < $tolerance)and(attr HARM < 0))
	vector do (QCOMP = BCOMP) (store5)
{}        
        vector do (segid = $sena) (store5)
        evaluate ($domain = $domain + 1)

  end loop newdomain

end if
display ------------------------------------------------------


! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! "best" rotational fit algorithms starts here
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

coordinates
	fit	
	selection = (segid AP1)
	mass = true
	lsq = true
end

! calculate COM and mass-based weigth

vector show sum (mass) (segid AP1)
evaluate ($MASS = $RESULT)
segment 
   	name = "DPT"
   	molecule name=DUPT number=1 end 
end
coordinates initialize
		selection = (segid DPT)
end	
vector do (X = 0) (segid DPT)
vector do (Y = 0) (segid DPT)
vector do (Z = 0) (segid DPT)
coordinates translate
	selection = (segid DPT)
	vector = (head = (segid AP1))
end
vector show element (X) (segid DPT)
evaluate ($XCOM = $RESULT)
vector show element (Y) (segid DPT)
evaluate ($YCOM = $RESULT)
vector show element (Z) (segid DPT)
evaluate ($ZCOM = $RESULT)

! info on AP0 structure (nonconverged rest)

vector show sum (mass) (segid AP0)
evaluate ($MASS9 = $RESULT)

! display some output

display
display ======================================================
display segid AP0 (rest)
display ======================================================
display MASS ..................	$MASS9
display >>>>> Residues:
for $1 in id (tag and (segid AP0AP0)) loop showloop
	vector show element (resid) (id $1)
	evaluate ($rnumber = $RESULT)
	vector show element (resname) (id $1)
	evaluate ($rname = $RESULT)
	display $rnumber $rname
end loop showloop	
display
display ======================================================
display segid AP1
display ======================================================
display MASS ..................	$MASS
display COM ...................	( $XCOM $YCOM $ZCOM )
display >>>>> Residues:
for $1 in id (tag and (segid AP1)) loop showloop
	vector show element (resid) (id $1)
	evaluate ($rnumber = $RESULT)
	vector show element (resname) (id $1)
	evaluate ($rname = $RESULT)
	display $rnumber $rname
end loop showloop	

evaluate ($domain = 1)
while ($domain < $ndomains) loop alldomains

	evaluate ($domain = $domain + 1)
	coordinates
		fit	
		selection = (segid AP1)
		mass = true
		lsq = true
	end
	evaluate ($sena = "AP" + encode ($domain))
	
	! calculate COM 
	
	vector show sum (mass) (segid $sena)
	evaluate ($MASS1 = $RESULT)
	coordinates swap
		selection = (segid $sena)
	end
	vector do (X = 0) (segid DPT)
	vector do (Y = 0) (segid DPT)
	vector do (Z = 0) (segid DPT)
	coordinates translate
		selection = (segid DPT)
		vector = (head = (segid $sena))
	end
	vector show element (X) (segid DPT)
	evaluate ($XCOM2 = $RESULT)
	vector show element (Y) (segid DPT)
	evaluate ($YCOM2 = $RESULT)
	vector show element (Z) (segid DPT)
	evaluate ($ZCOM2 = $RESULT)
	coordinates swap
		selection = (segid $sena)
	end
	vector do (X = 0) (segid DPT)
	vector do (Y = 0) (segid DPT)
	vector do (Z = 0) (segid DPT)
	coordinates translate
		selection = (segid DPT)
		vector = (head = (segid $sena))
	end
	vector show element (X) (segid DPT)
	evaluate ($XCOM1 = $RESULT)
	vector show element (Y) (segid DPT)
	evaluate ($YCOM1 = $RESULT)
	vector show element (Z) (segid DPT)
	evaluate ($ZCOM1 = $RESULT)
	
	! display some output
	
	display
	display ======================================================
	display segid $sena
	display ======================================================
	display MASS ..................	$MASS1
	display COM 1st ...............	( $XCOM1 $YCOM1 $ZCOM1 )
	display COM 2nd	...............	( $XCOM2 $YCOM2 $ZCOM2 )

	! calculate normalvector of bisecting plane
	
	evaluate ($XDD1 = $XCOM2 - $XCOM1)
	evaluate ($YDD1 = $YCOM2 - $YCOM1)
	evaluate ($ZDD1 = $ZCOM2 - $ZCOM1)
	evaluate ($CDIS = 	
		(sqrt($XDD1*$XDD1+$YDD1*$YDD1+$ZDD1*$ZDD1)))
	evaluate ($XPL = $XDD1 / $CDIS)
	evaluate ($YPL = $YDD1 / $CDIS)
	evaluate ($ZPL = $ZDD1 / $CDIS)

	! calculate bisecting point
	
	evaluate ($XBI = 0.5 * ($XCOM1 + $XCOM2))
	evaluate ($YBI = 0.5 * ($YCOM1 + $YCOM2))
	evaluate ($ZBI = 0.5 * ($ZCOM1 + $ZCOM2))

	! calculate best rmsd to compare

	coordinate fit
		mass = true
		selection = (segid $sena)
	end
	evaluate ($TH1 = $THETA1)
	evaluate ($TH2 = $THETA2)
	evaluate ($TH3 = $THETA3)
	coordinate rms 
		mass = true
		selection = (segid $sena)
	end
	evaluate ($RMSID = $RESULT)
	coordinate	
		fit	
		selection = (segid AP1)
		mass = true
		lsq = true
	end

	! calculate 2nd and 3rd point of plane (COM1 after rot about COM2)
	
	vector do (X = $XCOM1) (segid DPT)
	vector do (Y = $YCOM1) (segid DPT)
	vector do (Z = $ZCOM1) (segid DPT)
	coordinate rotate
		selection = (segid DPT)
		center = ( $XCOM2 $YCOM2 $ZCOM2 )
		euler = ( $TH1 $TH2 $TH3 )	
	end
	vector show element (X) (segid DPT)
	evaluate ($XP1 = $RESULT)
	vector show element (Y) (segid DPT)
	evaluate ($YP1 = $RESULT)
	vector show element (Z) (segid DPT)
	evaluate ($ZP1 = $RESULT)
	coordinate rotate
		selection = (segid DPT)
		center = ( $XCOM2 $YCOM2 $ZCOM2 )
		euler = ( $TH1 $TH2 $TH3 )	
	end
	vector show element (X) (segid DPT)
	evaluate ($XP2 = $RESULT)
	vector show element (Y) (segid DPT)
	evaluate ($YP2 = $RESULT)
	vector show element (Z) (segid DPT)
	evaluate ($ZP2 = $RESULT)

	! calculate axis of rotation

	evaluate ($XDEL1 = $XCOM1 - $XP1)
	evaluate ($YDEL1 = $YCOM1 - $YP1)
	evaluate ($ZDEL1 = $ZCOM1 - $ZP1)
	evaluate ($XDEL2 = $XCOM1 - $XP2)
	evaluate ($YDEL2 = $YCOM1 - $YP2)
	evaluate ($ZDEL2 = $ZCOM1 - $ZP2)
	evaluate ($XRIDEAL= ($YDEL2*$ZDEL1-$ZDEL2*$YDEL1))
	evaluate ($YRIDEAL= ($ZDEL2*$XDEL1-$XDEL2*$ZDEL1))
	evaluate ($ZRIDEAL= ($XDEL2*$YDEL1-$YDEL2*$XDEL1))
	evaluate ($NORG = 
		(sqrt($XRIDEAL*$XRIDEAL+$YRIDEAL*$YRIDEAL+$ZRIDEAL*$ZRIDEAL)))
	evaluate ($XRIDEAL= $XRIDEAL / $NORG)
	evaluate ($YRIDEAL= $YRIDEAL / $NORG)
	evaluate ($ZRIDEAL= $ZRIDEAL / $NORG)
	
	! calculate projection of COM2 on 3-point-plane

	evaluate ($PRID = ($XRIDEAL*$XDD1+$YRIDEAL*$YDD1+$ZRIDEAL*$ZDD1))
	evaluate ($XNEW = $XCOM2 - $PRID * $XRIDEAL)
	evaluate ($YNEW = $YCOM2 - $PRID * $YRIDEAL)
	evaluate ($ZNEW = $ZCOM2 - $PRID * $ZRIDEAL)
	
	! calculate least square angle
	
	evaluate ($XDTA1 = $XNEW - $XCOM1)
	evaluate ($YDTA1 = $YNEW - $YCOM1)
	evaluate ($ZDTA1 = $ZNEW - $ZCOM1)
	evaluate ($XDTA2 = $XNEW - $XP1)
	evaluate ($YDTA2 = $YNEW - $YP1)
	evaluate ($ZDTA2 = $ZNEW - $ZP1)
	evaluate ($NRG1 = 	
		(sqrt($XDTA1*$XDTA1+$YDTA1*$YDTA1+$ZDTA1*$ZDTA1)))
	evaluate ($NRG2 = 	
		(sqrt($XDTA2*$XDTA2+$YDTA2*$YDTA2+$ZDTA2*$ZDTA2)))
	evaluate ($COSINE = (($XDTA1 * $XDTA2 + $YDTA1 * $YDTA2 + 
		$ZDTA1 * $ZDTA2) / ($NRG1 * $NRG2)))
	evaluate ($ANGL = acos($COSINE))
	
	! calculate projection of rot axis on bisecting plane

	evaluate ($PERP = ($XRIDEAL*$XPL+$YRIDEAL*$YPL+$ZRIDEAL*$ZPL))
	evaluate ($ANGP = abs(asin($PERP)))	! angle of RIDEAL vs. PRO
	evaluate ($XPRO = $XRIDEAL - $PERP * $XPL)
	evaluate ($YPRO = $YRIDEAL - $PERP * $YPL)
	evaluate ($ZPRO = $ZRIDEAL - $PERP * $ZPL)
	evaluate ($NRM = 
		(sqrt($XPRO*$XPRO+$YPRO*$YPRO+$ZPRO*$ZPRO)))
	evaluate ($XPRO= $XPRO / $NRM)
	evaluate ($YPRO= $YPRO / $NRM)
	evaluate ($ZPRO= $ZPRO / $NRM)

	! calculate projected angle

	evaluate ($ANGLE = $ANGL*($XRIDEAL*$XPRO+$YRIDEAL*$YPRO+$ZRIDEAL*$ZPRO))
	
	! calculate hingepoint

	evaluate ($XDIR= ($YPL*$ZPRO-$ZPL*$YPRO))		
	evaluate ($YDIR= ($ZPL*$XPRO-$XPL*$ZPRO))
	evaluate ($ZDIR= ($XPL*$YPRO-$YPL*$XPRO))
	evaluate ($DIS = 0.5 * $CDIS / tan (0.5 * $ANGLE))
	evaluate ($XHI = $XBI + $DIS * $XDIR)
	evaluate ($YHI = $YBI + $DIS * $YDIR)
	evaluate ($ZHI = $ZBI + $DIS * $ZDIR)
		
	! calculate rmsd

	coordinate rotate
		axis = ( $XPRO $YPRO $ZPRO ) $ANGLE
		center = ( $XHI $YHI $ZHI )
	end
	coordinate rms 
		mass = true
		selection = (segid $sena)
	end
	evaluate ($RMSPRO = $RESULT)
	evaluate ($RELERR = (($RMSPRO / $RMSID) - 1)*100)
	coordinate	
		fit	
		selection = (segid AP1)
		mass = true
		lsq = true
	end

	! display some output
	display >>>>> Results:
	display ANGLE (hinge rot.) ....	$ANGLE deg.
	display HINGE ................. ( $XHI $YHI $ZHI )
	display ROT AXIS ..............	( $XPRO $YPRO $ZPRO )
	display >>>>> Accuracy:
	display COM-distance ..........	$CDIS A 	
	display Proj. angle ...........	$ANGP deg. (should be small for best results)
	display RMSD (least sq.) ......	$RMSID A
	display RMSD (proj)............	$RMSPRO A    
	display Relative error ........ $RELERR % 
	display >>>>> Residues:
	for $1 in id (tag and (segid $sena)) loop showloop
		vector show element (resid) (id $1)
		evaluate ($rnumber = $RESULT)
		vector show element (resname) (id $1)
		evaluate ($rname = $RESULT)
		display $rnumber $rname
	end loop showloop	

	! create dummy molecule to display vectors, axis
	
	evaluate ($duname = "DUM" + encode($domain))
	segment 
   		name = $duname
   		molecule name=DUM number=1 end 
	end
	coordinates initialize
		selection = (segid $duname)
	end	
	evaluate ($XAX1 = $XHI + 50 * $XPRO)
	evaluate ($YAX1 = $YHI + 50 * $YPRO)
	evaluate ($ZAX1 = $ZHI + 50 * $ZPRO)
	evaluate ($XAX2 = $XHI - 50 * $XPRO)
	evaluate ($YAX2 = $YHI - 50 * $YPRO)
	evaluate ($ZAX2 = $ZHI - 50 * $ZPRO)
	evaluate ($XAX3 = $XHI - 48 * $XPRO)
	evaluate ($YAX3 = $YHI - 48 * $YPRO)
	evaluate ($ZAX3 = $ZHI - 48 * $ZPRO)
	evaluate ($XDC1 = $XCOM1 - $XBI)
	evaluate ($YDC1 = $YCOM1 - $YBI)
	evaluate ($ZDC1 = $ZCOM1 - $ZBI)
	evaluate ($NRMC = 
		(sqrt($XDC1*$XDC1+$YDC1*$YDC1+$ZDC1*$ZDC1)))
	evaluate ($XDC1 = 2 * $XDC1 / $NRMC)
	evaluate ($YDC1 = 2 * $YDC1 / $NRMC)
	evaluate ($ZDC1 = 2 * $ZDC1 / $NRMC)
	evaluate ($XDC2 = $XCOM2 - $XBI)
	evaluate ($YDC2 = $YCOM2 - $YBI)
	evaluate ($ZDC2 = $ZCOM2 - $ZBI)
	evaluate ($NRMC = 
		(sqrt($XDC2*$XDC2+$YDC2*$YDC2+$ZDC2*$ZDC2)))
	evaluate ($XDC2 = 2 * $XDC2 / $NRMC)
	evaluate ($YDC2 = 2 * $YDC2 / $NRMC)
	evaluate ($ZDC2 = 2 * $ZDC2 / $NRMC)
	evaluate ($XPT1 = $XAX3 + $XDC1)
	evaluate ($YPT1 = $YAX3 + $YDC1)
	evaluate ($ZPT1 = $ZAX3 + $ZDC1)
	evaluate ($XPT2 = $XAX3 + $XDC2)
	evaluate ($YPT2 = $YAX3 + $YDC2)
	evaluate ($ZPT2 = $ZAX3 + $ZDC2)
	vector do (X = $XHI) ((segid $duname)and(name DU0))
	vector do (Y = $YHI) ((segid $duname)and(name DU0))
	vector do (Z = $ZHI) ((segid $duname)and(name DU0))
	vector do (X = $XCOM1) ((segid $duname)and(name DU1))
	vector do (Y = $YCOM1) ((segid $duname)and(name DU1))
	vector do (Z = $ZCOM1) ((segid $duname)and(name DU1))
	vector do (X = $XCOM2) ((segid $duname)and(name DU2))
	vector do (Y = $YCOM2) ((segid $duname)and(name DU2))
	vector do (Z = $ZCOM2) ((segid $duname)and(name DU2))
	vector do (X = $XAX1) ((segid $duname)and(name DU3))
	vector do (Y = $YAX1) ((segid $duname)and(name DU3))
	vector do (Z = $ZAX1) ((segid $duname)and(name DU3))
	vector do (X = $XAX2) ((segid $duname)and(name DU4))
	vector do (Y = $YAX2) ((segid $duname)and(name DU4))
	vector do (Z = $ZAX2) ((segid $duname)and(name DU4))
	vector do (X = $XAX3) ((segid $duname)and(name DU7))
	vector do (Y = $YAX3) ((segid $duname)and(name DU7))
	vector do (Z = $ZAX3) ((segid $duname)and(name DU7))
	vector do (X = $XPT1) ((segid $duname)and(name DU5))
	vector do (Y = $YPT1) ((segid $duname)and(name DU5))
	vector do (Z = $ZPT1) ((segid $duname)and(name DU5))
	vector do (X = $XPT2) ((segid $duname)and(name DU6))
	vector do (Y = $YPT2) ((segid $duname)and(name DU6))
	vector do (Z = $ZPT2) ((segid $duname)and(name DU6))

end loop alldomains

delete
	selection = (segid DPT)
end

write coordinates 
	output = $oname
end
write structure 
	output = $uname
end

stop