GNU Scientific Library、パッケージの多次元最小化を使用して、関数の最小値を見つけようとしています。私が使用している方法は、関数に実装されているBroyden-Fletcher-Goldfarb-Shanno(BFGS)アルゴリズムです。
gsl_multimin_fdfminimizer_vector_bfgs2
残念ながら、最初の反復後
gsl_multimin_fdfminimizer_iterate
コード27でエラーが発生しました:反復が解決に向けて進んでいません。
いくつかの初期点と異なる許容誤差/ステップの組み合わせを試しましたが、それでも同じエラーが表示されます。ここでの犯人を教えてください。
アイデアは、13個の原子のクラスターのポテンシャルエネルギーの最小値を見つけることです。エネルギーは関数my_fで計算されます。ここで、X、Y、Zは原子座標です。使用されるパラメータは
#define uA 0.1221
#define uXi 1.316
#define uP 8.612
#define uQ 2.516
#define uR0 2.8637
私が最小化しようとしている関数は、Gupta-アルミニウム13原子クラスターの可能性です。
//ここに私のグプタの可能性があります//
double
my_f (const gsl_vector *v, void *params)
{
double *p = (double *)params;
int ii,jj;
double vv = 0;
double sum1 = 0;
double sum2 = 0;
for(ii=1;ii<=NN;ii++){
X[ii]=gsl_vector_get (v, 3*ii-3);
Y[ii]=gsl_vector_get (v, 3*ii-2);
Z[ii]=gsl_vector_get (v, 3*ii-1 );
}
for(ii = 1; ii <= NN; ii++){
sum1 = 0;
for(jj = 1; jj <= NN; jj++){
if(jj!=ii){
sum1 += uA*exp(0 - uP*((rr(ii,jj)/uR0) - 1.));
}
}
sum2 = 0;
for(jj = 1; jj <= NN; jj++){
if(jj!=ii){
sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,jj)/uR0) - 1.));
}
}
vv += sum1 - sqrt(sum2);
}
return vv;
}
//ここに私のグラデーション関数があります//
double dVdr(int ii, int jj){
double vr;
double sum1 = 0;
double sum2 = 0;
int kk;
vr = 0. - 2.*uA*(uP/uR0)*exp(0. - uP*((rr(ii,jj)/uR0) - 1.));
for(kk = 1; kk <= NN; kk++){
if(kk != ii){
sum1 += uXi*uXi*exp(0 - 2*uQ*((rr(ii,kk)/uR0) - 1.));
}
if(kk != jj){
sum2 += uXi*uXi*exp(0 - 2*uQ*((rr(kk,jj)/uR0) - 1.));
}
}
vr -= 0.5*uXi*uXi*(0. - 2.*uQ/uR0)*exp(0. - 2.*uQ*((rr(ii,jj)/uR0) - 1.))*( (1./sqrt(sum1)) + (1./sqrt(sum2)) );
return vr;
}
//-----------------------------------------------
void Force(int jj, double & fx, double & fy, double & fz){
int mm;
double dxx, dyy, dzz, r;
fx = 0;
fy = 0;
fz = 0;
for(mm = 1; mm <= NN; mm++){
if(mm != jj){
dxx = X[mm] - X[jj];
dyy = Y[mm] - Y[jj];
dzz = Z[mm] - Z[jj];
r = sqrt(dxx*dxx + dyy*dyy + dzz*dzz);
fx += dVdr(mm,jj)*(dxx)/r;
fy += dVdr(mm,jj)*(dyy)/r;
fz += dVdr(mm,jj)*(dzz)/r;
}
}
}
//--------------------------------------------
void
my_df (const gsl_vector *v, void *params,
gsl_vector *df)
{
int ii;
double fx,fy,fz;
double *p = (double *)params;
for(ii=1;ii<=NN;ii++){
X[ii]=gsl_vector_get (v, 3*ii-3);
Y[ii]=gsl_vector_get (v, 3*ii-2);
Z[ii]=gsl_vector_get (v, 3*ii-1 );
}
for(ii=1;ii<=NN;ii++){
Force(ii,fx,fy,fz);
gsl_vector_set(df, 3*ii-3, fx );
gsl_vector_set(df, 3*ii-2, fy );
gsl_vector_set(df, 3*ii-1, fz );
}
}
//----------------------------------