Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

437 lines (367 sloc) 17.082 kb
/***********************************************************************
*
* Copyright (C) 2008 Thomas Chiarappa, Carsten Urbach
*
* This file is part of tmLQCD.
*
* tmLQCD is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* tmLQCD is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with tmLQCD. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#ifdef HAVE_CONFIG_H
# include<config.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "global.h"
#include "su3.h"
#include "su3adj.h"
#include "linalg_eo.h"
#include "start.h"
#include "linsolve.h"
#include "solver/solver.h"
#include "deriv_Sb.h"
#include "tm_operators.h"
#include "chebyshev_polynomial.h"
#include "Nondegenerate_Matrix.h"
#include "Hopping_Matrix.h"
#include "phmc.h"
#include "Nondegenerate_Matrix.h"
#include "chebyshev_polynomial_nd.h"
#include "Ptilde_nd.h"
#include "reweighting_factor_nd.h"
#include "monomial.h"
#include "hamiltonian_field.h"
#include "ndpoly_monomial.h"
extern int phmc_exact_poly;
/********************************************
*
* Here \delta S_b is computed
*
********************************************/
void ndpoly_derivative(const int id, hamiltonian_field_t * const hf) {
int j, k;
monomial * mnl = &monomial_list[id];
/* This factor 2 a missing factor 2 in trace_lambda */
(*mnl).forcefactor = -2.*phmc_Cpol*phmc_invmaxev;
/* Recall: The GAMMA_5 left of delta M_eo is done in deriv_Sb !!! */
if (g_epsbar!=0.0 || phmc_exact_poly==0){
/* Here comes the definitions for the chi_j fields */
/* from j=0 (chi_0 = phi) ..... to j = n-1 */
/* in g_chi_up_spinor_field[0] (g_chi_dn_spinor_field[0] we expect */
/* to find the phi field, the pseudo fermion field */
/* i.e. must be equal to mnl->pf (mnl->pf2) */
assign(g_chi_up_spinor_field[0], mnl->pf, VOLUME/2);
assign(g_chi_dn_spinor_field[0], mnl->pf2, VOLUME/2);
for(k = 1; k < (phmc_dop_n_cheby-1); k++) {
Q_tau1_min_cconst_ND(g_chi_up_spinor_field[k], g_chi_dn_spinor_field[k],
g_chi_up_spinor_field[k-1], g_chi_dn_spinor_field[k-1],
phmc_root[k-1]);
}
/* Here comes the remaining fields chi_k ; k=n,...,2n-1 */
/*They are evaluated step-by-step overwriting the same field (phmc_dop_n_cheby)*/
assign(g_chi_up_spinor_field[phmc_dop_n_cheby], g_chi_up_spinor_field[phmc_dop_n_cheby-2], VOLUME/2);
assign(g_chi_dn_spinor_field[phmc_dop_n_cheby], g_chi_dn_spinor_field[phmc_dop_n_cheby-2], VOLUME/2);
for(j=(phmc_dop_n_cheby-1); j>=1; j--) {
assign(g_chi_up_spinor_field[phmc_dop_n_cheby-1], g_chi_up_spinor_field[phmc_dop_n_cheby], VOLUME/2);
assign(g_chi_dn_spinor_field[phmc_dop_n_cheby-1], g_chi_dn_spinor_field[phmc_dop_n_cheby], VOLUME/2);
Q_tau1_min_cconst_ND(g_chi_up_spinor_field[phmc_dop_n_cheby], g_chi_dn_spinor_field[phmc_dop_n_cheby],
g_chi_up_spinor_field[phmc_dop_n_cheby-1], g_chi_dn_spinor_field[phmc_dop_n_cheby-1],
phmc_root[2*phmc_dop_n_cheby-3-j]);
/* Get the even parts of the (j-1)th chi_spinors */
H_eo_ND(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1], EO);
/* \delta M_eo sandwitched by chi[j-1]_e^\dagger and chi[2N-j]_o */
deriv_Sb(EO, g_spinor_field[DUM_DERI], g_chi_up_spinor_field[phmc_dop_n_cheby], hf); /* UP */
deriv_Sb(EO, g_spinor_field[DUM_DERI+1], g_chi_dn_spinor_field[phmc_dop_n_cheby], hf); /* DN */
/* Get the even parts of the (2N-j)-th chi_spinors */
H_eo_ND(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
g_chi_up_spinor_field[phmc_dop_n_cheby], g_chi_dn_spinor_field[phmc_dop_n_cheby], EO);
/* \delta M_oe sandwitched by chi[j-1]_o^\dagger and chi[2N-j]_e */
deriv_Sb(OE, g_chi_up_spinor_field[j-1], g_spinor_field[DUM_DERI], hf);
deriv_Sb(OE, g_chi_dn_spinor_field[j-1], g_spinor_field[DUM_DERI+1], hf);
}
}
else if(g_epsbar == 0.0) {
/* Here comes the definitions for the chi_j fields */
/* from j=0 (chi_0 = phi) ..... to j = n-1 */
assign(g_chi_up_spinor_field[0], mnl->pf, VOLUME/2);
for(k = 1; k < (phmc_dop_n_cheby-1); k++) {
Qtm_pm_min_cconst_nrm(g_chi_up_spinor_field[k],
g_chi_up_spinor_field[k-1],
phmc_root[k-1]);
}
assign(g_chi_up_spinor_field[phmc_dop_n_cheby],
g_chi_up_spinor_field[phmc_dop_n_cheby-2], VOLUME/2);
for(j = (phmc_dop_n_cheby-1); j >= 1; j--) {
assign(g_chi_up_spinor_field[phmc_dop_n_cheby-1],
g_chi_up_spinor_field[phmc_dop_n_cheby], VOLUME/2);
Qtm_pm_min_cconst_nrm(g_chi_up_spinor_field[phmc_dop_n_cheby],
g_chi_up_spinor_field[phmc_dop_n_cheby-1],
phmc_root[2*phmc_dop_n_cheby-3-j]);
Qtm_minus_psi(g_spinor_field[DUM_DERI+3],g_chi_up_spinor_field[j-1]);
H_eo_tm_inv_psi(g_spinor_field[DUM_DERI+2], g_chi_up_spinor_field[phmc_dop_n_cheby], EO, -1.);
deriv_Sb(OE, g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI+2], hf);
H_eo_tm_inv_psi(g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+3], EO, 1.);
deriv_Sb(EO, g_spinor_field[DUM_DERI+2], g_chi_up_spinor_field[phmc_dop_n_cheby], hf);
Qtm_minus_psi(g_spinor_field[DUM_DERI+3],g_chi_up_spinor_field[phmc_dop_n_cheby]);
H_eo_tm_inv_psi(g_spinor_field[DUM_DERI+2],g_spinor_field[DUM_DERI+3], EO, +1.);
deriv_Sb(OE, g_chi_up_spinor_field[j-1] , g_spinor_field[DUM_DERI+2], hf);
H_eo_tm_inv_psi(g_spinor_field[DUM_DERI+2], g_chi_up_spinor_field[j-1], EO, -1.);
deriv_Sb(EO, g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI+3], hf);
}
}
/*
Normalisation by the largest EW is done in update_momenta
using mnl->forcefactor
*/
}
void ndpoly_heatbath(const int id, hamiltonian_field_t * const hf) {
int j;
double temp;
monomial * mnl = &monomial_list[id];
(*mnl).energy0 = 0.;
random_spinor_field(g_chi_up_spinor_field[0], VOLUME/2, (*mnl).rngrepro);
(*mnl).energy0 = square_norm(g_chi_up_spinor_field[0], VOLUME/2, 1);
if(g_epsbar!=0.0 || phmc_exact_poly == 0){
random_spinor_field(g_chi_dn_spinor_field[0], VOLUME/2, (*mnl).rngrepro);
(*mnl).energy0 += square_norm(g_chi_dn_spinor_field[0], VOLUME/2, 1);
}
else {
zero_spinor_field(g_chi_dn_spinor_field[0], VOLUME/2);
}
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Here comes the computation of H_old with \n \n");
printf("PHMC: First: random spinors and their norm \n ");
printf("PHMC: OLD Ennergy UP %e \n", (*mnl).energy0);
printf("PHMC: OLD Energy DN + UP %e \n\n", (*mnl).energy0);
}
if(phmc_exact_poly==0){
QNon_degenerate(g_chi_up_spinor_field[1], g_chi_dn_spinor_field[1],
g_chi_up_spinor_field[0], g_chi_dn_spinor_field[0]);
for(j = 1; j < (phmc_dop_n_cheby); j++){
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
assign(g_chi_dn_spinor_field[0], g_chi_dn_spinor_field[1], VOLUME/2);
Q_tau1_min_cconst_ND(g_chi_up_spinor_field[1], g_chi_dn_spinor_field[1],
g_chi_up_spinor_field[0], g_chi_dn_spinor_field[0],
phmc_root[phmc_dop_n_cheby-2+j]);
}
Poly_tilde_ND(g_chi_up_spinor_field[0], g_chi_dn_spinor_field[0], phmc_ptilde_cheby_coef,
phmc_ptilde_n_cheby, g_chi_up_spinor_field[1], g_chi_dn_spinor_field[1]);
}
else if( phmc_exact_poly==1 && g_epsbar!=0.0) {
/* Attention this is Q * tau1, up/dn are exchanged in the input spinor */
/* this is used as an preconditioner */
QNon_degenerate(g_chi_up_spinor_field[1],g_chi_dn_spinor_field[1],
g_chi_dn_spinor_field[0],g_chi_up_spinor_field[0]);
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
assign(g_chi_dn_spinor_field[0], g_chi_dn_spinor_field[1], VOLUME/2);
/* solve Q*tau1*P(Q^2) *x=y */
cg_her_nd(g_chi_up_spinor_field[1],g_chi_dn_spinor_field[1],
g_chi_up_spinor_field[0],g_chi_dn_spinor_field[0],
1000,1.e-16,0,VOLUME/2, Qtau1_P_ND);
/* phi= Bdagger phi */
for(j = 1; j < (phmc_dop_n_cheby); j++){
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
assign(g_chi_dn_spinor_field[0], g_chi_dn_spinor_field[1], VOLUME/2);
Q_tau1_min_cconst_ND(g_chi_up_spinor_field[1], g_chi_dn_spinor_field[1],
g_chi_up_spinor_field[0], g_chi_dn_spinor_field[0],
phmc_root[phmc_dop_n_cheby-2+j]);
}
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
assign(g_chi_dn_spinor_field[0], g_chi_dn_spinor_field[1], VOLUME/2);
}
else if(phmc_exact_poly==1 && g_epsbar==0.0) {
Qtm_pm_psi(g_chi_up_spinor_field[1], g_chi_up_spinor_field[0]);
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
/* solve (Q+)*(Q-)*P((Q+)*(Q-)) *x=y */
cg_her(g_chi_up_spinor_field[1], g_chi_up_spinor_field[0],
1000,1.e-16,0,VOLUME/2, Qtm_pm_Ptm_pm_psi);
/* phi= Bdagger phi */
for(j = 1; j < (phmc_dop_n_cheby); j++){
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
Qtm_pm_min_cconst_nrm(g_chi_up_spinor_field[1],
g_chi_up_spinor_field[0],
phmc_root[phmc_dop_n_cheby-2+j]);
}
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
}
assign(mnl->pf, g_chi_up_spinor_field[0], VOLUME/2);
assign(mnl->pf2, g_chi_dn_spinor_field[0], VOLUME/2);
temp = square_norm(g_chi_up_spinor_field[0], VOLUME/2, 1);
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Then: evaluate Norm of pseudofermion heatbath BHB \n ");
printf("PHMC: Norm of BHB up squared %e \n", temp);
}
if(g_epsbar!=0.0 || phmc_exact_poly==0)
temp += square_norm(g_chi_dn_spinor_field[0], VOLUME/2, 1);
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)){
printf("PHMC: Norm of BHB up + BHB dn squared %e \n\n", temp);
}
if(g_proc_id == 0 && g_debug_level > 3) {
printf("called ndpoly_heatbath for id %d with g_running_phmc = %d\n", id, g_running_phmc);
}
return;
}
double ndpoly_acc(const int id, hamiltonian_field_t * const hf) {
int j, ij=0;
double temp, sgn, fact, Diff;
double Ener[8];
double factor[8];
monomial * mnl = &monomial_list[id];
spinor *up0, *dn0, *up1, *dn1, *dummy;
mnl->energy1 = 0.;
Ener[0] = 0;
factor[0] = 1.0;
for(j = 1; j < 8; j++){
factor[j] = j*factor[j-1];
Ener[j] = 0;
}
/* IF PHMC */
up0 = g_chi_up_spinor_field[0];
up1 = g_chi_up_spinor_field[1];
dn0 = g_chi_dn_spinor_field[0];
dn1 = g_chi_dn_spinor_field[1];
/* This is needed if we consider only "1" in eq. 9 */
assign(up0, mnl->pf , VOLUME/2);
assign(dn0, mnl->pf2, VOLUME/2);
if(phmc_exact_poly==0) {
for(j = 1; j <= (phmc_dop_n_cheby-1); j++) {
/* Change this name !!*/
Q_tau1_min_cconst_ND(up1, dn1, up0, dn0, phmc_root[j-1]);
dummy = up1; up1 = up0; up0 = dummy;
dummy = dn1; dn1 = dn0; dn0 = dummy;
/* result always in up0 and dn0 */
}
ij=0;
if(up0 != g_chi_up_spinor_field[ij]) {
assign(g_chi_up_spinor_field[ij], up0, VOLUME/2);
assign(g_chi_dn_spinor_field[ij], dn0, VOLUME/2);
}
temp = square_norm(g_chi_up_spinor_field[ij], VOLUME/2, 1);
Ener[ij] = temp;
temp = square_norm(g_chi_dn_spinor_field[ij], VOLUME/2, 1);
Ener[ij] += temp;
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Here comes the computation of H_new with \n \n");
printf("PHMC: At j=%d PHMC Final Energy %e \n", ij, mnl->energy1+Ener[ij]);
printf("PHMC: At j=%d PHMC Only Final Energy %e \n", ij, Ener[ij]);
}
/* Here comes the loop for the evaluation of A, A^2, ... */
for(j = 1; j < 1; j++){ /* To omit corrections just set j<1 */
if(j % 2){ /* Chi[j] = ( Qdag P Ptilde ) Chi[j-1] */
Poly_tilde_ND(g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j],
phmc_ptilde_cheby_coef, phmc_ptilde_n_cheby,
g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1]);
QdaggerQ_poly(g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1],
phmc_dop_cheby_coef, phmc_dop_n_cheby,
g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j]);
QdaggerNon_degenerate(g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j],
g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1]);
}
else { /* Chi[j] = ( Ptilde P Q ) Chi[j-1] */
QNon_degenerate(g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j],
g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1]);
QdaggerQ_poly(g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1],
phmc_dop_cheby_coef, phmc_dop_n_cheby, g_chi_up_spinor_field[j],
g_chi_dn_spinor_field[j]);
Poly_tilde_ND(g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j],
phmc_ptilde_cheby_coef, phmc_ptilde_n_cheby,
g_chi_up_spinor_field[j-1], g_chi_dn_spinor_field[j-1]);
}
Ener[j] = Ener[j-1] + Ener[0];
sgn = -1.0;
for(ij = 1; ij < j; ij++){
fact = factor[j] / (factor[ij] * factor[j-ij]);
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Here j=%d and ij=%d sign=%f fact=%f \n", j ,ij, sgn, fact);
}
Ener[j] += sgn*fact*Ener[ij];
sgn = -sgn;
}
temp = square_norm(g_chi_up_spinor_field[j], VOLUME/2, 1);
temp += square_norm(g_chi_dn_spinor_field[j], VOLUME/2, 1);
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Here j=%d sign=%f temp=%e \n", j, sgn, temp);
}
Ener[j] += sgn*temp;
Diff = fabs(Ener[j] - Ener[j-1]);
if((g_proc_id == g_stdio_proc) && (g_debug_level > 0)) {
printf("PHMC: Correction aftern %d steps: %e \n", j, Diff);
}
if(Diff < g_acc_Hfin) {
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: At j = %d PHMC Only Final Energy %e \n", j, Ener[j]);
}
break;
}
}
mnl->energy1 += Ener[ij]; /* this is quite sticky */
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: At j = %d P=%e +HMC Final Energy %e \n\n", ij, Ener[ij], mnl->energy1);
}
}
else if(phmc_exact_poly==1 && g_epsbar!=0.0) {
/* B(Q*tau1) */
for(j = 1; j <= (phmc_dop_n_cheby-1); j++){
Q_tau1_min_cconst_ND(up1, dn1, up0, dn0, phmc_root[j-1]);
dummy = up1; up1 = up0; up0 = dummy;
dummy = dn1; dn1 = dn0; dn0 = dummy;
/* result always in up0 and dn0 */
}
if(up0 != g_chi_up_spinor_field[0]) {
assign(g_chi_up_spinor_field[0], up0, VOLUME/2);
assign(g_chi_dn_spinor_field[0], dn0, VOLUME/2);
}
temp = square_norm(g_chi_up_spinor_field[0], VOLUME/2, 1);
Ener[0] = temp;
temp = square_norm(g_chi_dn_spinor_field[0], VOLUME/2, 1);
Ener[0] += temp;
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
ij=0;
printf("PHMC: Here comes the computation of H_new with \n \n");
printf("PHMC: At j=%d P+HMC Final Energy %e \n", ij, mnl->energy1+Ener[0]);
printf("PHMC: At j=%d PHMC Only Final Energy %e \n", ij, Ener[0]);
}
mnl->energy1 += Ener[0];
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: At j = %d P=%e +HMC Final Energy %e \n\n", ij, Ener[0], mnl->energy1);
}
}
else if(phmc_exact_poly == 1 && g_epsbar == 0.0) {
for(j = 1; j < (phmc_dop_n_cheby); j++) {
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
Qtm_pm_min_cconst_nrm(g_chi_up_spinor_field[1],
g_chi_up_spinor_field[0],
phmc_root[j-1]);
}
assign(g_chi_up_spinor_field[0], g_chi_up_spinor_field[1], VOLUME/2);
temp = square_norm(g_chi_up_spinor_field[0], VOLUME/2, 1);
Ener[0] = temp;
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: Here comes the computation of H_new with \n \n");
printf("PHMC: At j=%d P+HMC Final Energy %e \n", ij, mnl->energy1+Ener[0]);
printf("PHMC: At j=%d PHMC Only Final Energy %e \n", ij, Ener[0]);
}
mnl->energy1 += Ener[0];
if((g_proc_id == g_stdio_proc) && (g_debug_level > 2)) {
printf("PHMC: At j = %d P=%e +HMC Final Energy %e \n\n", ij, Ener[0], mnl->energy1);
}
}
if(g_proc_id == 0 && g_debug_level > 3) {
printf("called ndpoly_acc for id %d %d dH = %1.4e\n", id, g_running_phmc, mnl->energy1 - mnl->energy0);
}
/* END IF PHMC */
return(mnl->energy1 - mnl->energy0);
}
Jump to Line
Something went wrong with that request. Please try again.