//   LMOD minimization

#include "xmin_opt.h"
#include "lmod_opt.h"

// LMOD minmization:

struct xmin_opt xo;
struct lmod_opt lo;


molecule mol;
int natm;
int eig_recalc, ndim_arnoldi;
int lmod_restart, n_best_conf,mc_option;
int n_lmod_steps, nlig, apply_rigdock;
int n_pose_try, random_seed, verbosity, error_flag;
float ene, minim_grms, ene_window, rtemp;
int lig_start[ dynamic ], lig_end[ dynamic ], lig_cent[ dynamic ];
float x[ dynamic ], g[ dynamic ], fret,
      conflib[ dynamic ], lmod_traj[ dynamic ];
float tr_min[ dynamic ], tr_max[ dynamic ], rot_min[ dynamic ],
      rot_max[ dynamic ];
float dgrad;
//float surften;
int i,max_pdb_files;
float glob_min_energy;
float grms;
string fname,sys_command;
point dummy;
int ier;


    mol = getpdb( "gb1_c2_ab_lmod_1.pdb" );
    readparm( mol, "c2_ab.prmtop" );
    natm = mol.natoms;

    allocate x[ 3*natm ]; allocate g[ 3*natm ];

    setxyz_from_mol( mol, NULL, x );

    xmin_opt_init( xo );

    xo.maxiter   = 100000;
    xo.grms_tol  = 0.05;
    xo.method    = 3;
    xo.m_lbfgs   = 3;
    xo.print_level = 2;

//    dielc=1.0;
//    mm_options( "ntpr=999, gb=1, rgbmax=999.0, cut=999.0, scee=1.0, scnb=1.0, diel=C");
//    mm_options( "ntpr=1, gb=1, gbsa=1, surften=0.0072, kappa=0.10395, cut=999.0, diel=C ");
    mm_options( "ntpr=10000, gb=1, gbsa=1, surften=0.0072, kappa=0.0, rgbmax=25.0, cut=999.0, scee=1.2, scnb=1.2, diel=C, epsext=78.3, temp0=298.15, e_debug=2");
//    mm_options( "ntpr=999, gb=1, rgbmax=999.0, cut=999.0, scee=1.0, scnb=1.0, diel=C, dielc=4.0 ");
    mme_init( mol, NULL, "::ZZZ", dummy, NULL );

    ene = mme( x, g, 0 );
    ene = xmin( mme, natm, x, g, ene, grms, xo);
        
    if( mytaskid == 0 ) {
  
      printf( " Unfrozen minimized energy = %12.3lf\n" , ene );
      printf( " Graduebt RMS              = %12.3lf\n" , grms );
      printf( " Iter                      = %12d\n"    , xo.iter );
      printf( " CPU time                  = %12.3f\n\n", xo.xmin_time );
    }

// Newton-Raphson minimization
//    mm_options( "ntpr=1" );
    
    dgrad = 3*natm*0.001;
    newton( x, 3*natm, fret, mme, mme2, dgrad, 0.0, 15 );

    setmol_from_xyz( mol, NULL, x );                            // load x[] into mol

    fname = sprintf( "gb1_ab_nmode.pdb" );
    putpdb( fname, mol, "-brook" );                             // save conf in pdb file

for(i=1; i <= 3*natm; i=i+6)
{
    printf("%12.7f%12.7f%12.7f%12.7f%12.7f%12.7f\n", x[i], x[i+1], x[i+2], x[i+3], x[i+4], x[i+5]);
}

// get the normal modes
    nmode( x, 3*natm, mme2, 0, 0, 0.0, 0.0, 0);