function [Vtot_opt,Vtot_pJ] = V(q)
% computes V for MonteCarlo problem.  the first output is used in the
% minimization of the potential and is V in pJ times 1e12.  The second is
% what is used in the MC steps and is return in picoJoules.

% O1 = (0,0,0)
% O2 = (xO2,0,0)
% H1 = (xH1,yH1,0);

xO2=q(1);  %all x,y,z in Angstroms
xH1=q(2); yH1=q(3);
xH2=q(4); yH2=q(5); zH2=q(6);

% calculate the variables used in the potentials.  Included is converting
% cartesian coordinate into angles
ROO=abs(xO2);
ROH1=norm([xH1; yH1]);
ROH2=norm([xH2-xO2; yH2; zH2]);

Rlong1=norm([xH1-xO2; yH1]);
HOO=acos(-0.5*(Rlong1*Rlong1-(ROH1*ROH1+ROO*ROO))./(ROH1*ROO));

Rlong2=norm([xH2; yH2; zH2]);
OOH=acos(-0.5*(Rlong2*Rlong2-(ROH2*ROH2+ROO*ROO))./(ROH2*ROO));

cosphi=yH1*yH2/(abs(yH1)*norm([yH2; zH2]));

kOO=3e-6;  % [=]picoJoule/Angstrom^2
LO=1.6;  % [=] Angstrom

% stretching potential... also utilizes the Morse potential for OH bonds
Vstretch=Vmorse(ROH1)+Vmorse(ROH2)+0.5*kOO*(ROO-LO).^2;

% angular potential function
ktheta=1e-6;%[=]picoJoule/radian^2
theta0=1.8; %[=]radian
Vbend=0.5*ktheta*(((HOO-theta0).^2)+((OOH-theta0).^2));

% dihedral potential function
Kphi=80e3/6.022e11; % picoJoule (per molecule)
Vphi=Kphi*sin(0.5*HOO)*sin(0.5*OOH)*(cosphi-cos(1.7)).^2;

% Total potential is the sum of all three
Vtot_pJ = Vstretch + Vbend + Vphi;

Vtot_opt = Vtot_pJ*1e12;

%q0 =
%    1.6000
%   -0.2386
%    1.0225
%    1.8386
%   -0.1317
%    1.0140
%V(q0)
%
%ans =
%
%  1.8982e-034