Cholesky  decomposition  method  for  canonical  molecular  orbital  calculations  of  proteinsToshiyuki  Hirano  and  Fumitoshi  Sato  1I n t r o d u c t i o n   Institute  of  Industrial  Science,  the  University  of  Tokyo,  Japan    The   canonical   molecular   orbital   calculations   of   proteins  provide   significant   theoretical   views,   which   are   difficult   to  observe   experimentally,   and   are   promising   technology   in  various   fields   such   as   drug   discovery   and   nanoelectronics  devices  [1].  Even  with  modern  high-­‐performance  computers  and   parallel   technologies   in   quantum   chemistry,   however,  electronic   structure   calculations   of   proteins   still   remain   a  time-­‐consuming   and   laborious   task.   We   aimed   at   the  development   of   computational   methods   to   efficiently  calculate   electronic   structure   of   proteins   in   distributed  parallel   computing,   which   becomes   mainstream   now.   In  particular,  efficiency  of  two-­‐electron  repulsion  integral  (ERI)  computations   based   on   Cholesky   decomposition   (CD)  method  is  focused  in  this  study.  
The  rest  of  paper  is  constructed  as  follows.  In  the  section  2,  we   discuss   the   theoretical   foundation   based   on   the   CD  method   and   the   use   of   the   pivoted   Cholesky   decomposition  to   calculate   low-­‐rank   approximations   of   matrices.   The  section  3  is  devoted  to  numerical  results  using  ProteinDF  on  distributed   memory   parallel   computer.   Finally,   we   discuss  capabilities   related   to   the   analyses   of   electronic   structure  calculations  of  proteins  in  the  section  4.  
2

M e t h o d o l o g y   2.1In  self-­‐consistent  field  (SCF)  calculations,  the  construction  of  the  Fock  matrix  is  one  of  the  most  time-­‐consuming  steps.  In  particular,   much   computer   time   is   spent   on   the   task   of  computing  the  Coulomb  part  𝑱  and  the  Fock  exchange  part  𝑲,  which   include   ERI   calculations.   The   matrix   𝑱   and   𝑲   are  formally  computed  as  follows:  𝑃!" 𝑝𝑞 𝑟𝑠   (1)  𝑃!! 𝑝𝑞 𝑟𝑠   (2)  !"𝐾!" =!"

Here,  𝑷  denotes   the   density   matrix,   and   the   ERI   𝑝𝑞 𝑟𝑠   defined  in  Eq.  (3).  
1(3)  𝑝𝑞 𝑟𝑠 =𝑔!∗ 𝒓 𝑔! 𝒓
𝑔∗ 𝒓! 𝑔! 𝒓! d𝒓d𝒓!   
𝒓 − 𝒓! !
The   super   matrix  𝑉!",!" ,   which   is   defined   by  𝑉!",!" = 𝑝𝑞 𝑟𝑠 ,  is   symmetric   and   positive   semi-­‐definite,   and   can   be  decomposed   by   means   of   the   Cholesky   procedure   into   a  product  of  a  lower  triangular  matrix  𝑳  as  follows:  !

𝑉!",!" = 𝑝𝑞 𝑟𝑠 ≈𝐿!,!" 𝐿!,!"   !!!

(4)  The   accuracy   of   the   Cholesky   representation   of   the   ERI,   Eq.  
(4),  is  measured  by  the  residual.  Consequently,  the  accuracy  can  be  controlled  by  decomposition  threshold,  𝛿.    

(5)  𝐿!,!" 𝐿!,!" ≤ 𝛿  !!!

According   to   some   computational   results,   the   Cholesky  decomposition   (CD)   method   gives   reliable,   and   balanced  results   even   with   rather   high   decomposition   thresholds  (𝛿 = 1.0×10!! )  [2].    
The   number   of   dimensions   of   the   original   super   matrix  𝑉  formally   becomes   𝑁 𝑁 + 1 /2   with   𝑁   basis   functions,   so  such   huge   memory   for   the   matrix   may   not   be   assigned   in  each   computer   node.   Therefore,   the   size   of   the   matrix   is  reduced   in   advance   of   the   CDAM   method   [3].   The   diagonal  elements  of  the  𝑽  matrix,  𝑉!",!" = 𝑝𝑞 𝑝𝑞 ,  are  preliminarily  estimated  and  screened  by  threshold  𝜏.  
𝑝𝑞 𝑝𝑞 ≈ 𝐼 ! 𝐼 ! ≥ 𝜏  (6)  Here,  the  AO  products,  𝐼 = 𝑝𝑞,  can  be  regarded  as  candidates  of   the   Cholesky   basis.   Note   that   the   screening   technique   of  the  Eq.  (6)  is  associated  with  the  Cauchy-­‐Schwarz  screening,  which  is  very  useful  for  providing  rigorous  ERI  bound.  Thus,  the   number   of   non-­‐negligible   basis   function   pairs   will   scale  as  𝒪 𝑁 ,  thanks  to  locality  of  basis  functions.    Once  the  Cholesky  vectors  𝑳  are  obtained,  the  Coulomb  part  𝑱  and  exchange  part  𝑲  of  the  Fock  matrix  are  approximately  gained  as:  𝐽!" =Cholesky  decomposition  method  𝐽!" =!

𝑉!",!" −(𝑝𝑞, 𝑟𝑠)𝑃!" ≈!"𝐾!" =12!"

𝐿!,!" 𝐿!,!" 𝑃!"   (7)  𝑋!,!" 𝑋!,!"   (8)  !

(𝑝𝑟, 𝑞𝑠)𝑃!" ≈!"𝑋!,!" =!

!

𝐿!,!" 𝑄!"   (9)  ! 𝑄!,! 𝑄!,!     (10)  !

𝑃!" =In  this  study,  the  matrix  𝑿  is  generated  by  𝑸,  which  is  given  by  the  CD  of  the  density  matrix  𝑷,  instead  of  MO  coefficient  matrix.  Then,  we  can  use  the  density  matrix  as  initial  guess  for  the  SCF  calculation.  
Once  𝑳  is   estimated,   the   procedures   in   the  𝑱  and  𝑲  matrices  generation  can  be  given  by  simple  matrix-­‐matrix  operations  without  the  ERI  estimations.  Since  the  well-­‐optimized  linear  algebra  library,  such  as  BLAS,  for  various  computer  systems  is   available,   the  𝑱  and  𝑲  matrices   are   constructed   by   simple  and  effective  method.  
2.2Low-­‐rank  approximate  decomposition  Cholesky  To   get   smaller   super   matrix   without   losing   calculation  precision   is   an   advantage   in   terms   of   memory   capacity   and  number  of  floating-­‐point  arithmetic  operations.  The  CD  of  n-­‐dimensional   matrix   usually   gives   n-­‐dimensional   triangular  matrix.   In   this   study,   the   pivoted   CD   to   compute   low-­‐rank  approximations   [4]   is   adopted.   The   resulting   truncation  error  is  rigorously  controlled  in  terms  of  the  trace  norm.    In  general,  accessing  all  elements  of  the  full  matrix  is  needed  in   CD.   Consequently,   the   huge   full   matrix   should   be  temporarily   stored   in   memory.   The   elements   that   are   left  unused  in  the  pivoted  CD  are  also  included  in  the  full  matrix.  
In  order  to  produce  Cholesky  vectors  on  limited  memory,  it  is  desired  that  the  full  matrix  is  not  kept  in  memory,  and  that  the   ERIs,   which   are   elements   of   the   full   matrix,   are   directly  calculated.   Putting   all   the   above   components   together,   we  arrive  at  the  algorithm  (Fig.  1).  Since  there  is  no  necessary  to  obtain   the   huge   super   matrix   temporarily,   the   Cholesky  vectors  are  gained  by  minimal  computer  resources.  
Although  the  Cholesky  vectors  are  directly  obtained  and  the  huge   super   matrix   is   not   required,   the   Cholesky   vectors  remains  large  size,  and  may  not  be  stored  in  memory.  Hence  the   Cholesky   vectors   𝑳   should   be   distributed   to   all  processesor  elements  (PEs).  In  this  study,  the  entries  of  the  Cholesky   vectors   are   stored   as   set   of   row   vectors   for  efficiency.   On   this   algorithm,   the   access   direction   of   the  matrix   is   only   the   row   direction.   With   keeping   diagonal  elements,  ℓ𝓁!,!! ,   in   all   PEs,   the   operations   for   row-­‐vectors  elements   of   the   Cholesky   vectors   held   in   each   PE   are  independent  (line:  11),  so  no  communications  are  required.  
And,  the  molecular  integrals  may  be  carried  out  only  for  the  row-­‐vectors  of  Cholesky  vectors  (line:  9)。This  method  can  be  directly  parallel  processing  the  Cholesky  vectors.  
3R e s u l t s    a n d    d i s c u s s i o n   We  implemented  this  method  to  ProteinDF  program  [5].  The  canonical   molecular   orbital   calculations   of   Insulin   (51  residues)  and  Interleukin  (133  residues)  are  carried  out  by  using   ProteinDF.   The   number   of   AO   pairs   in   Insulin   and  Interleukin,   which   were   9,952,491   and   70,918,095,   were  dramatically  reduced  to  169,842(1.7%)  and  460,785  (0.6%)  by   CDAM   pre-­‐screening   (𝜏 = 10!! ),   respectively.   So,   the  combination  of  CDAM  and  pivoted  CD  were  considered  to  be  valid   and   feasible   for   canonical   molecular   orbital  computations.   In   order   to   clarify   effective   threshold   value,  further  investigation  is  required.  
4the   method   gives   smaller   and   optimized   Cholesky   vectors.  
We  emphasis  that  the  Cholesky  vectors  𝑳  are  calculated  and  stored  before  the  SCF  loop,  the  Coulomb  and  exchange  part  of   Fock   matrix   are   given   by   the   multiplication   of   matrices,  without   expensive   ERI   computations.   Though   required  memory  is  larger  than  the  direct  method,  in  which  the  ERIs  are   calculated   at   every   SCF   iteration   in   return,   we   can   get  high   efficiency   of   computation   on   the   distributed   parallel  computer   systems.   The   memory   size   per   computer   node  tends   to   increase,   and   this   problem   will   be   overcome.   We  suppose   that   this   method   is   useful   for   electronic   structure  calculation   of   the   large   molecule   system   in   parallel   and  distributed  memory  system.    
A c k n o w l e d g e m e n t s   A   part   of   this   research   was   done   in   “Revolutionary  Simulation   Software   for   the   21st   century   (RSS21)”   project  supported   by   Research   and   Development   for   Next-­‐
generation  Information  Technology  of  Ministry  of  Education,  Culture,   Sports,   Science   and   Technology   (MEXT).   And   this  research   was   partly   supported   by   “Research   and  Development   of   the   Next-­‐Generation   Integrated   Simulation  of  Living  Matter”,  a  part  of  the  Development  and  Use  of  the  Next-­‐Generation   Supercomputer   Project   of   the   Ministry   of  Education,  Culture,  Sports,  Science  and  Technology  (MEXT).  
We   thank   the   Institute   of   Statistical   Mathematics   for   the  facilities  and  the  use  of  Fujitsu  PRIMERGY  RX200S5.  
R e f e r e n c e s   [1]

[2]

[3]

[4]

C o n c l u s i o n   A   new   ERI   parallel   computation   method   based   on   modified  CDAM   method   and   the   pivoted   low-­‐approximate   CD   was  developed  and  implemented  to  ProteinDF.  The  method  was  able   to   control   the   efficiency   of   molecular   integral  estimations   with   improving   the   accuracy.   As   the   way   of  achievement   on   the   ERI   computation   efficiency   in   this  method  resembles  the  well-­‐known  cutoff  method,  negligible  ERI   values   are   strictly   estimated   based   on   CD.   Additionally,  [5]Y.   Tokita,   S.   Yamada,   W.   Luo,   Y.   Goto,   N.   Bouley-­‐Ford,   H.  
Nakajima,   Y.   Watanabe,   Protein   Photoconductors   and  Photodiodes,  Angew.  Chem.  Int.  Ed.,  50,  11663  –11666,  2011.  
F.   Aquilante,   R.   Lindh,   T.B.   Pedersen,   Unbiased   auxiliary   basis  sets  for  accurate  two-­‐electron  integral  approximations,  J.  Chem.  
Phys.,  127,  114107,  2007.  
Y.  Okiyam,  T.  Nakano,  K.  Yamashita,  Y.  Mochizuki,  N.  Taguchi,  S.  
Tanaka,   Acceleration   of   fragment   molecular   orbital  calculations  with  Cholesky  decomposition  approach,  Chemical  Physics  Letters,  490,  84-­‐89,  2010.  
H.   Harbrechta,   M.   Petersa,   R.   Schneider,   On   the   low-­‐rank  approximation   by   the   pivoted   Cholesky   decomposition,   Appl.  
Numer.  Math.,  62,  428-­‐440,  2012.  
T.   Inaba,   F.   Sato,   Development   of   parallel   density   functional  program   using   distributed   matrix   to   calculate   all-­‐electron  canonical   wavefunction   of   large   molecules,   J.   Comput.   Chem.,  28,  984-­‐995,  2007.  

line  1  2  3  4  5  6  pseudo  code  begin    set  𝑚 ∶=   0;    set  𝒅 ∶=   𝑑𝑖𝑎𝑔 𝑨   and  𝑒𝑟𝑟𝑜𝑟 =    𝒅 ! ;    initialize  𝝅 ∶=    1, 2, … , 𝑛 ;    while  𝑒𝑟𝑟𝑜𝑟 >   𝜀  do      set  𝑖 ∶=   𝑎𝑟𝑔𝑚𝑎𝑥 𝑑!! : 𝑗 = 𝑚, 𝑚 + 1, … , 𝑛 ;  7  8  9  10  11  12  13  14  15                  end                  swap  𝜋!   and  𝜋! ;  set  ℓ𝓁!,!! ≔ 𝑑!! ;  compute  molecular  integral  𝑎!! ,!! = 𝜋! 𝜋!    𝑖; 𝑚 + 1: 𝑛   for  𝑚 + 1 ≤ 𝑖 ≤ 𝑛  do    compute  ℓ𝓁!,!! ≔ 𝑎!! ,!! − !!!
!!! ℓ𝓁!,!! ℓ𝓁!,!! /ℓ𝓁!,!! ;    update  𝑑!! ≔ 𝑑!! − ℓ𝓁!!,!! ;  compute  𝑒𝑟𝑟𝑜𝑟 ≔ !!!!!! 𝑑!! ;  increase  𝑚 ≔ 𝑚 + 1;  Fig.  1.  Pseudo  code  for  pivoted  Cholesky  decomposition  of  ERIs