cellopt/CrystalSystems_t.cpp at master · tgruene/cellopt · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

#include "CrystalSystems.h"
#include <cmath>

/**
 * tetragonal cell: a = b and alpha = beta = gamma = 90; vary a, c
 * @param cell
 * @param restraints
 * @return
 */
double CrystSyst::Tetragonal::target(const gsl_vector* cell, void* restraints) {
     gsl_vector* R = (gsl_vector*) restraints;

        const double a (gsl_vector_get(cell, 0));
        const double c (gsl_vector_get(cell, 1));
    double T(0.0);
    for (size_t i = 0; i < R->size; i+=5) {
        const double dx (gsl_vector_get(R, i+0));
        const double dy (gsl_vector_get(R, i+1));
        const double dz (gsl_vector_get(R, i+2));
        const double R2 (gsl_vector_get(R, i+3));
        const double wR (gsl_vector_get(R, i+4));
        double dX2 = (dx*dx+dy*dy)*a*a + dz*dz*c*c;
        T += wR*(dX2 - R2)*(dX2 - R2);
    }
    return T;
}

void CrystSyst::Tetragonal::gradT(const gsl_vector* cell, void* restraints, gsl_vector* gradT) {
    gsl_vector* R = (gsl_vector*) restraints;
    const double a(gsl_vector_get(cell, 0));
    const double c(gsl_vector_get(cell, 1));

    double dTda(0), dTdc(0);
    for (size_t i = 0; i < R->size; i += 5) {
        const double dx(gsl_vector_get(R, i + 0));
        const double dy(gsl_vector_get(R, i + 1));
        const double dz(gsl_vector_get(R, i + 2));
        const double R2(gsl_vector_get(R, i + 3));
        const double wR(gsl_vector_get(R, i + 4));

        double dX2 = (dx*dx+dy*dy)*a*a + dz*dz*c*c;

        double dta = 2 * (dx * dx + dy * dy)* a;
        double dtc = 2 * dz * dz * c;

        dTda += 2 * wR * (dX2 - R2) * dta;
        dTdc += 2 * wR * (dX2 - R2) * dtc;
    }
    gsl_vector_set(gradT, 0, dTda);
    gsl_vector_set(gradT, 1, dTdc);

}

void CrystSyst::Tetragonal::TgradT(const gsl_vector* cell, void* restraints,
double* T, gsl_vector* gradT) {
    gsl_vector* R = (gsl_vector*) restraints;

    const double a (gsl_vector_get(cell, 0));
    const double c (gsl_vector_get(cell, 1));

    *T = 0.0;
    double dTda(0), dTdc(0);
    for (size_t i = 0; i < R->size; i += 5) {
        const double dx(gsl_vector_get(R, i + 0));
        const double dy(gsl_vector_get(R, i + 1));
        const double dz(gsl_vector_get(R, i + 2));
        const double R2(gsl_vector_get(R, i + 3));
        const double wR(gsl_vector_get(R, i + 4));

        double dX2 = (dx * dx + dy * dy)* a * a + dz * dz * c * c;

        double dta = 2 * (dx * dx + dy * dy)* a;
        double dtc = 2 * dz * dz * c;

        const double dX2R2 (dX2 - R2);
        *T       += wR*dX2R2*dX2R2;
        dTda     += 2 * wR*dX2R2 * dta;
        dTdc     += 2 * wR*dX2R2 * dtc;
    }
    gsl_vector_set(gradT, 0, dTda);
    gsl_vector_set(gradT, 1, dTdc);
}