Mam funkcję readPDB, która ma wczytać plik. Funkcja przyjmuje 2 argumenty - nazwa pliku oraz tablica, do której przypisywane są dane z pliku. Program jednak nie chcę się skompilować i wyświetla się błąd jak w temacie.
#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <ctime>
#include <cmath>
#include <iomanip>
#include <cstring>
#define N 10
using namespace std;
void writePDB(double r[][3], char nazwaPDB[])
{
int i;
fstream plik_pdb;
plik_pdb.open("AR_MD.pdb", ios::out | ios::app);
if(plik_pdb == NULL)
{
perror ("NIE MOGE OTWORZYC PLIKU PDB! ABY DOPISAC DO NIEGO DANE");
}
for(i=0;i<N;i++)
plik_pdb << "ATOM" << " " << fixed << setprecision(4) << i+1 << " " << "AR" << " " << "AR" << fixed << setprecision(7) << i+1 << " "
<< fixed << setprecision(3) << r[i][0] << "\t" << fixed << setprecision(3) << r[i][1] << "\t" << fixed << setprecision(3) << r[i][2] << "\t"
<< " 1.00 0.00 AAA AR" << endl;
plik_pdb << "END" << endl;
plik_pdb.close();
}
void readPDB(string nazwapliku, double tab[][3])
{
fstream plik;
plik.open(nazwapliku.c_str() );
int licznik=-1, w=0;
char napis[15];
if(plik == NULL){
perror ("NIE MOGE OTWORZYC PLIKU, BY POBRAC DANE!");
}
while (plik.good() == true){
plik >> napis;
if((napis[0]=='A')&&(napis[1]=='T')&&(napis[2]=='O')&&(napis[3]=='M')) licznik=0;
switch(licznik){
case 6:
tab[w][0]= atof(napis);
break;
case 7:
tab[w][1]= atof(napis);
break;
case 8:
tab[w][2]= atof(napis);
w++;
break;
default:
break;
}
licznik++;
}
plik.close();
}
double Ekin(double v[N][3],float Mas) //funkcja do mierzenia sredniej Energii kin
{
int i;
double Ek=0; //energia kinetyczna
for(i=0;i<N;i++)
Ek+=Mas*0.5*(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
return Ek;
}
double Eelec(double r[N][3], double qi[N], double qj[N]) //energia elektrostatyczna dla ladunkow losowych
{
int i,j;
double pi=M_PI, epsilon0=8.85*pow(10.0,-19.0), Eelec;
for(i=1;i<N;i++)
for(j=1;j<N;j++)
{
qi[i]=rand()%21-10;
qj[j]=rand()%21-10;
Eelec=qi[i]*qj[j]/4*pi*epsilon0*r[i][j];
}
return Eelec;
}
void printVel(double v[N][3])
{
int i;
fstream plikV;
plikV.open("Vel.dat", ios::out | ios::app);
if(plikV == NULL)
{
perror ("NIE MOGE OTWORZYC PLIKU PDB! ABY DOPISAC DO NIEGO DANE");
}
for(i=0;i<N;i++)
plikV << v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2] << endl;
plikV.close();
}
void randVel(double v[N][3], double vMax)
{
int i, j;
for (i = 0; i < N; i++)
for (j = 0; j < 3; j++)
v[i][j] = (vMax * 2.0 * rand() / (RAND_MAX - 1.0))-vMax; //predkooci na razie losowo
}
void Acc(const double r[N][3],double a[N][3], float Mas, float sigma, float eps)
{
int i,j;
double F;
double dx, dy, dz, dr;
for (i = 0; i < N; i++)
for (j = 0; j < 3; j++)
a[i][j] = 0; //zerowanie sil
for (i = 0; i < N; i++)
{
for (j = i +1; j < N; j++)
{
dx=0;
dy=0;
dz=0;
dx = r[i][0] - r[j][0];
dy = r[i][1] - r[j][1];
dz = r[i][2] - r[j][2];
dr=sqrt(dx*dx+dy*dy+dz*dz);
F= eps*(48 * pow(dr/sigma, -13.0) - 24 * pow(dr/sigma, -7.0));
a[i][0]+=F*dx/(Mas * dr);
a[j][0]-=F*dx/(Mas * dr);
a[i][1]+=F*dy/(Mas * dr);
a[j][1]-=F*dy/(Mas * dr);
a[i][2]+=F*dz/(Mas * dr);
a[j][2]-=F*dz/(Mas * dr);
}
}
}
int main()
{
double r[N][3]; // polozenia N-indeksuje atomy, x->0,y->1,z->2
double v[N][3]; // predkosci N-indeksuje atomy, x->0,y->1,z->2
double a[N][3]; // przyspieszenia N-indeksuje atomy, x->0,y->1,z->2
int i, j, k;
double vMax=5.775; // A/ps
double dt=0.001; // timestep krok czasowy w ps
int TIME=10000; //-number of timesteps
float L=14; //PBC rozmiar skrzynki 2*L x 2*L x 2*L
float sigma=3.81638, eps=0.23846; //parametryVDW
float Mas= 39.95; //MASA
char nPDB[20]="AR_MD.pdb"; //nazwa pliku PDB
double kb=1.38064852; //Stala boltzmanna w J/K
double Ek; //energia kinetyczna
double Ec; //energia elektrostatyczna
double qi[N], qj[N];
fstream plik_pdb, data, Ecsave;
srand(time(NULL));
cout << "Dynamika molekularna Argonu na podstawie algorytmu velocity Verlet rozszerzona o obliczenie oddzialywan elektrostatycznych " << endl;
cout << endl;
cout << "W pliku Multi Frame PDB zapisywane sa trajektorie wynikowe" << endl;
cout << endl;
cout << "W pliku data z rozszerzeniem dat zapisywane sa informacje o kroku czasowym, energii kinetycznej" << endl;
cout << endl;
cout << "W pliku Ecsave z rozszerzeniem dat zapisywane sa informacje o kroku czasowym i oddzialywaniach elektrostatycznych" << endl;
plik_pdb.open("AR_MD.pdb", ios::out | ios::app);
if(plik_pdb == NULL)
{
perror("NIE MOGE OTWORZYC PLIKU PDB!");
return 1;
}
plik_pdb.close();
data.open("data.dat", ios::out | ios::app);
if(plik_pdb == NULL)
{
perror("NIE MOGE OTWORZYC PLIKU PDB!");
return 1;
}
readPDB(nPDB, r[N][3]);
randVel(v,vMax); //generowanie losowych predkosci
writePDB(r,nPDB); //zapis pierwszych koordynatow
Ecsave.open("Ec.dat", ios::out | ios::app);
for(k=1;k<=TIME;k++){
Acc(r,a, Mas, sigma, eps);
for (i = 0; i < N; i++)
for (j = 0; j < 3; j++)
{
r[i][j] += v[i][j] * dt + 0.5 * a[i][j] * dt * dt;
// PBC
if (r[i][j] < -L) r[i][j] += 2*L;
if (r[i][j] > L) r[i][j] -= 2*L;
v[i][j] += 0.5 * a[i][j] * dt;
}
Acc(r,a, Mas, sigma, eps);
for (i = 0; i < N; i++)
for (j = 0; j < 3; j++)
v[i][j] += 0.5 * a[i][j] * dt;
if(k%10==0)
{
writePDB(r,nPDB);
Ek=Ekin(v,Mas);
Ec=Eelec(r, qi, qj);
//cout << Ec << endl;
data << k << fixed << setprecision(6) << "\t" << Ek << fixed << setprecision(6) << "\t" << 2.0*1.66054*Ek/(3*N*kb) << "\t"
<< endl;
Ecsave << k << "\t" << Ec << endl;
//2.0*1.66054*Ek/(3*N*kb) bierze sie z przeliczenia EK na J (1.66054) oraz podstawienia kb w J/K (1.38064852) 10^34 w liczniku i mianowniku skraca sie
}
}
Ecsave.close();
data.close();
cout << "GOTOWE!" << endl;
system("Pause");
return 0;
}