[DFTB-Plus-User] Lattice vectors initialization

[Riccardo Petraglia]

Hello everybody,

I am trying to implement an interface between I-Pi 
(http://arxiv.org/abs/1402.1045, https://github.com/i-pi/i-pi) and DFTB+ 
in order to use DFTB+ with enhanced MD methods.

The interface needs to add a way to insert atom positions and lattice 
vectors directly inside the code (without the needing of writing an 
input on a file) and reading the forces and stress tensor .

I am inserting the positions and lattice vectors just after the 
"lpGeomOpt" loop and reading the forces just before the first check on 
the result of the geometry optimization.

Everithing works fine if the lattice vectors inserted by my modification 
are the same (or really close: <2A) to the lattice vectors specified in 
the gen file. When the lattice vectors I insert are different from the 
lattice vectors specified in the gen file the SCC iteration never converge.

I understand it is kind of impossible debbugging a code by mail... That 
is not my purpose. I would like to know if the lattice vectors are used 
to initiliaze something before the "injection point" I chosen (taking 
into account that I am trying to spread the new values with the lines 
following the injection). Another way to ask the question is: when using 
the barostat, the change in the lattice vectors how are spread to every 
subroutine? Is it normal that above a certain tolerance value, about 2A 
per side in a cubic box, the SCC is not able to converge anymore?

Here below the lines I inserted in the dftb+.F90 file with some lines 
that should be helpful in determining the point to determine the 
position in relation to the original code. I used the version of DFTB+ 
1.2.1.
As attached the dftb_in.hsd and the geometry file.


   ipi = .true.
### Injecting the coordinates and lattice vectors

   lpGeomOpt: do while ((iGeoStep <= nGeoSteps) .or. ipi)

      write(81,*) Coord0, latvec #debug
      !ipi: Reading positions from ipi
      if (ipi) then
         call coordsFromSocket(natom, Coord0, latvec)         #sending 
the number of atom to the subroutine and asking for coordinate and 
latticevectors
         recVec2p = 
latVec(:,:)                                             #copied from 
some other part of the DFTB+ code in order to spread the latVec to all 
the subroutines
         call 
matinv(recVec2p)                                              #
         recVec2p = reshape(recVec2p, (/3, 3/), order=(/2, 
1/))             #
         recVec = 2.0_dp * pi * recVec2p          #
         volume = determinant33(latVec)          #
         recVolume = determinant33(recVec)           #
         if (tSCC) then          #
            print*, 'updating lattice vectors in 
SCC'                       #
            call updateLatVecs_SCC(latVec, recVec, 
volume)                  #
            mCutoff = max(mCutoff, 
getSCCCutoff())                          #
         end 
if                                                             #
         if (tDispersion) 
then                                              #
            print*, 'updating lattice vectors in dispersion 
modules'        #
            call updateLatVecs(myDispersion, latVec, recVec, volume, 
&      #
&specie0)                                                  #
            mCutoff = max(mCutoff, 
getRCutoff(myDispersion))                #
         end 
if                                                             #
      end if           #



### Reading forces and strees tensor

      if (ipi) then
         print*, 'Sending Forces to the socket'
         call forcesToSocket(natom, energy%ETotal, -1*totalDeriv, 
totalStress)
         cycle lpGeomOpt
      end if

       !! If geometry minimizer finished and the last calculated 
geometry is the
       !! minimal one (not necessary the case, depends on the optimizer!)
       !! -> we are finished.
       !! Otherwise we have to recalc everything in the converged geometry.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: c6.gen
Type: chemical/x-genbank
Size: 2780 bytes
Desc: not available
URL: <http://mailman.zfn.uni-bremen.de/pipermail/dftb-plus-user/attachments/20141112/4b3e0521/attachment.gen>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: c6.xyz
Type: chemical/x-xyz
Size: 2825 bytes
Desc: not available
URL: <http://mailman.zfn.uni-bremen.de/pipermail/dftb-plus-user/attachments/20141112/4b3e0521/attachment.xyz>
-------------- next part --------------
Geometry = GenFormat { 
     <<< "c6.gen"
} 
 
Driver = ConjugateGradient{ 
   ConvergentForcesOnly = No 
   MaxForceComponent = 1E-10
   MaxSteps = 50
   } 
Hamiltonian = DFTB{ 
   MaxSCCIterations = 100
   Charge = 0.0 
   DampXH = Yes 
   DampXHExponent = 4.0 
   Eigensolver = RelativelyRobust{} 
   HubbardDerivs = { 
      H = -0.1857 
      C = -0.1492
      S = -0.11
   } 
   MaxAngularMomentum = { 
      H = "s" 
      C = "p"
      S = "d"
   } 
   SCC = Yes 
   SlaterKosterFiles = Type2FileNames{ 
      Prefix = ../SK-parameters/3ob-1-1/ 
      Separator = "-" 
      Suffix = ".skf" 
   } 
   ThirdOrderFull = Yes 
   KPointsAndWeights[absolute] = {
      0.0 0.0 0.0 1.0
   }
} 
Options { 
   WriteResultsTag = Yes  
}