*/
/*
- There's something wrong with energy of the last residue
+ There's something wrong with energy calculations of redidues with E ≈ -0
*/
*
*/
+ if (CalculateAtomicDistances(
+ Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
+ >= minimalCAdistance)
+ {
+ return void();
+ }
+
const float kCouplingConstant = 27.888;
const float minimalAtomDistance{ 0.5 },
minEnergy{ -9.9 };
float HbondEnergy{ 0 };
float distanceNO{ 0 }, distanceHC{ 0 }, distanceHO{ 0 }, distanceNC{ 0 };
- std::cout << "For Donor №" << Donor.info->nr - 1 << " and accpetor №" << Acceptor.info->nr - 1 << std::endl;
-
- if( !(Donor.is_proline) ){
- if (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
- && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || (initParams.addHydrogens) ) )
- {
-
- distanceNO = CalculateAtomicDistances(
- Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
- distanceNC = CalculateAtomicDistances(
- Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
-
- if (initParams.addHydrogens){
- if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
- rvec atomH{};
- float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
- for (int i{XX}; i <= ZZ; ++i){
- float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
- atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
- atomH[i] += prevCO / prevCODist;
- }
- distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
- distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+ std::cout << "For Donor №" << Donor.info->nr - 1 << " and Accpetor №" << Acceptor.info->nr - 1 << std::endl;
+
+ if( !(Donor.is_proline) && (Acceptor.getIndex(backboneAtomTypes::AtomC) && Acceptor.getIndex(backboneAtomTypes::AtomO)
+ && Donor.getIndex(backboneAtomTypes::AtomN) && ( Donor.getIndex(backboneAtomTypes::AtomH) || initParams.addHydrogens ) ) ){ // TODO
+ distanceNO = CalculateAtomicDistances(
+ Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+ distanceNC = CalculateAtomicDistances(
+ Donor.getIndex(backboneAtomTypes::AtomN), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+ if (initParams.addHydrogens){
+ if (Donor.prevResi != nullptr && Donor.prevResi->getIndex(backboneAtomTypes::AtomC) && Donor.prevResi->getIndex(backboneAtomTypes::AtomO)){
+ rvec atomH{};
+ float prevCODist {CalculateAtomicDistances(Donor.prevResi->getIndex(backboneAtomTypes::AtomC), Donor.prevResi->getIndex(backboneAtomTypes::AtomO), fr, pbc)};
+ for (int i{XX}; i <= ZZ; ++i){
+ float prevCO = fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomC)][i] - fr.x[Donor.prevResi->getIndex(backboneAtomTypes::AtomO)][i];
+ atomH[i] = fr.x[Donor.getIndex(backboneAtomTypes::AtomH)][i]; // Но на самом деле берутся координаты N
+ atomH[i] += prevCO / prevCODist;
}
- else{
- distanceHO = distanceNO;
- distanceHC = distanceNC;
- }
- }
- else {
- distanceHO = CalculateAtomicDistances(
- Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
- distanceHC = CalculateAtomicDistances(
- Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
- }
-
- if (CalculateAtomicDistances(
- Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc)
- < minimalCAdistance)
- {
- if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
- || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
- {
- HbondEnergy = minEnergy;
- }
- else
- {
- HbondEnergy =
- kCouplingConstant
- * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
-// HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
-
-// std::cout.precision(5);
-// std::cout << "Calculated ENERGY = " << HbondEnergy << std::endl;
-
-// if ( HbondEnergy == 0){
-// std::cout << "Calculated ENERGY = " << HbondEnergy << " For donor " << Donor.info->nr << " and acceptor " << Acceptor.info->nr << std::endl;
-// }
-
- if ( HbondEnergy < minEnergy ){
- HbondEnergy = minEnergy;
- }
- }
- }
+ distanceHO = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+ distanceHC = CalculateAtomicDistances(atomH, Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+ }
+ else{
+ distanceHO = distanceNO;
+ distanceHC = distanceNC;
+ }
+ }
+ else {
+ distanceHO = CalculateAtomicDistances(
+ Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomO), fr, pbc);
+ distanceHC = CalculateAtomicDistances(
+ Donor.getIndex(backboneAtomTypes::AtomH), Acceptor.getIndex(backboneAtomTypes::AtomC), fr, pbc);
+ }
+ if ((distanceNO < minimalAtomDistance) || (distanceHC < minimalAtomDistance)
+ || (distanceHO < minimalAtomDistance) || (distanceNC < minimalAtomDistance))
+ {
+ HbondEnergy = minEnergy;
+ }
+ else{
+ HbondEnergy =
+ kCouplingConstant
+ * ((1 / distanceNO) + (1 / distanceHC) - (1 / distanceHO) - (1 / distanceNC));
}
-
-
std::cout << "CA-CA distance: " << CalculateAtomicDistances(
Donor.getIndex(backboneAtomTypes::AtomCA), Acceptor.getIndex(backboneAtomTypes::AtomCA), fr, pbc) << std::endl;
std::cout << "H-O distance: " << distanceHO << std::endl;
std::cout << "H-C distance: " << distanceHC << std::endl;
- }
-
+ HbondEnergy = std::round(HbondEnergy * 1000) / 1000;
- std::cout << "Calculated energy = " << HbondEnergy << std::endl;
+ if ( HbondEnergy < minEnergy ){
+ HbondEnergy = minEnergy;
+ }
+ std::cout << "Calculated energy = " << HbondEnergy << std::endl;
+ }
+ else{
+ std::cout << "Donor Is Proline" << std::endl;
+ }
- if (HbondEnergy < Donor.acceptorEnergy[0]){
- Donor.acceptor[1] = Donor.acceptor[0];
- Donor.acceptor[0] = Acceptor.info;
- Donor.acceptorEnergy[0] = HbondEnergy;
- }
- else if (HbondEnergy < Donor.acceptorEnergy[1]){
- Donor.acceptor[1] = Acceptor.info;
- Donor.acceptorEnergy[1] = HbondEnergy;
- }
+ if (HbondEnergy < Donor.acceptorEnergy[0]){
+ Donor.acceptor[1] = Donor.acceptor[0];
+ Donor.acceptor[0] = Acceptor.info;
+ Donor.acceptorEnergy[0] = HbondEnergy;
+ }
+ else if (HbondEnergy < Donor.acceptorEnergy[1]){
+ Donor.acceptor[1] = Acceptor.info;
+ Donor.acceptorEnergy[1] = HbondEnergy;
+ }
- if (HbondEnergy < Acceptor.donorEnergy[0]){
- Acceptor.donor[1] = Acceptor.donor[0];
- Acceptor.donor[0] = Donor.info;
- Acceptor.donorEnergy[0] = HbondEnergy;
- }
- else if (HbondEnergy < Acceptor.donorEnergy[1]){
- Acceptor.donor[1] = Donor.info;
- Acceptor.donorEnergy[1] = HbondEnergy;
- }
+ if (HbondEnergy < Acceptor.donorEnergy[0]){
+ Acceptor.donor[1] = Acceptor.donor[0];
+ Acceptor.donor[0] = Donor.info;
+ Acceptor.donorEnergy[0] = HbondEnergy;
+ }
+ else if (HbondEnergy < Acceptor.donorEnergy[1]){
+ Acceptor.donor[1] = Donor.info;
+ Acceptor.donorEnergy[1] = HbondEnergy;
+ }
}