#include<iostream>
#include "ATOMIC_3x9.h"
#include <stdio.h>
#include <iomanip>

using namespace std;

int main(void){

 ATOMIC VPOT;

 char direc[1024]; //char file[1024];

 //THESE ARE THE SLIGHTLY MODIFIED DFT DATA FOR THE ATOMIC PES
 //SUITABLE FOR CRP METHOD
 sprintf(direc,"/home/gernot/H2Pt211/Atomic/DFT_ATOMIC_MOD");

 //INITIALIZE THE PES
 VPOT.INIT_ATOMIC(direc);
 //CALL THE PES
 //VPOT.CALC_ATOMIC(pos_cart, dvd_cart, &vpot);


 ofstream out; out.open("test.pot");

 double LX = 6.955149254;
 double LY = 2.839427793;
 double ang_a0 = 1./0.52917721;
 double zmin = -3.65; double zmax = 7.05;

 double vpot = 0.;
 double pos_cart[3]; double dvd_cart[3];
 pos_cart[0] = LX/9. * 0. * ang_a0; 
 pos_cart[1] = LY/4. * 0. * ang_a0;


 for(int i = 0; i < 109; i++){
 
   pos_cart[2] = 0.1* i + zmin;
   pos_cart[2] *= ang_a0;   
   VPOT.CALC_ATOMIC(pos_cart, dvd_cart, &vpot);  
   out << setprecision(10) << showpoint << pos_cart[2]/ang_a0 << " " << (vpot)*27.211383-75.150547 << "\n";

 }

 //optimize Z to obatin minimum binding energies
  ofstream outi; outi.open("opt_pes.dat");
  for(int i = 0; i < 100; i++){
    pos_cart[0] = LX/99. * i * ang_a0;
    for(int j = 0; j < 100; j++){
      pos_cart[1] = LY/99. * j * ang_a0;
      double VREF = 1600.; double ZREF = 0.;
      double dzz = 4.5/1599. ;
      for(int k = 0; k < 1600; k++){
        pos_cart[2] = -1.5 + dzz*k * ang_a0;
        VPOT.CALC_ATOMIC(pos_cart, dvd_cart, &vpot);
        if(vpot < VREF){
         ZREF = pos_cart[2];
         VREF = vpot;
        }
      }
      outi << showpoint << setprecision(10) << pos_cart[0]/ang_a0 << " " << pos_cart[1]/ang_a0 << " " << ZREF/ang_a0 << " " << VREF << "\n";
    }
    outi << "\n";
  } 
  outi.close();
  


}