phosphorene nanotube optimization

[beauyy]

Dear sir,
     I built a phosphorene nanotube and optimized it. But during the optimization process, the structure was difficult to converge and all the P atoms were discreted from the bulk. T The .py file is attached as following. What's wrong with my optimization parameter?

# -------------------------------------------------------------                           
# Bulk configuration                                                                       
# -------------------------------------------------------------                           
                                                                                           
# Set up lattice                                                                           
vector_a = [44.2622, 0.0, 0.0]*Angstrom                                                   
vector_b = [0.0, 44.2622, 0.0]*Angstrom                                                   
vector_c = [0.0, 0.0, 3.31586338739]*Angstrom                                             
lattice = UnitCell(vector_a, vector_b, vector_c)                                           
                                                                                           
# Define elements                                                                         
elements = [Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus,                   
            Phosphorus, Phosphorus, Phosphorus, Phosphorus, Phosphorus]                   
                                                                                           
# Define coordinates                                                                       
fractional_coordinates = [[ 0.36582515,  0.48499453,  0.25000081],                         
                          [ 0.47283227,  0.3677228 ,  0.25000081],                         
                          [ 0.61743134,  0.43325366,  0.25000081],                         
                          [ 0.59979135,  0.5910257 ,  0.25000081],                         
                          [ 0.44429017,  0.62300331,  0.25000081],                         
                          [ 0.44429017,  0.37699669,  0.74999917],                         
                          [ 0.59979135,  0.4089743 ,  0.74999917],                         
                          [ 0.61743134,  0.56674634,  0.74999917],                         
                          [ 0.47283227,  0.6322772 ,  0.74999917],                         
                          [ 0.36582515,  0.51500547,  0.74999917],                         
                          [ 0.40783012,  0.34090592,  0.74999943],                         
                          [ 0.62284893,  0.36314607,  0.74999943],                         
                          [ 0.66814176,  0.57451369,  0.74999943],                         
                          [ 0.48111546,  0.68290592,  0.74999943],                         
                          [ 0.32023402,  0.53852838,  0.74999943],                         
                          [ 0.37712092,  0.36338569,  0.24999903],                         
                          [ 0.59197975,  0.34088652,  0.24999903],                         
                          [ 0.67977275,  0.53827678,  0.24999903],                         
                          [ 0.51917299,  0.68276983,  0.24999903],                         
                          [ 0.33212387,  0.57468119,  0.24999903],                         
                          [ 0.40044368,  0.40909823,  0.24999996],                         
                          [ 0.55571165,  0.37719373,  0.24999996],                         
                          [ 0.63403508,  0.51500332,  0.24999996],                         
                          [ 0.52717365,  0.63207883,  0.24999996],                         
                          [ 0.38280622,  0.56662589,  0.24999996],                         
                          [ 0.38280623,  0.43337411,  0.75000002],                         
                          [ 0.52717365,  0.36792117,  0.75000002],                         
                          [ 0.63403508,  0.48499668,  0.75000002],                         
                          [ 0.55571165,  0.62280626,  0.75000002],                         
                          [ 0.40044368,  0.59090176,  0.75000002],                         
                          [ 0.33212387,  0.42531881,  0.75000097],                         
                          [ 0.51917299,  0.31723017,  0.75000097],                         
                          [ 0.67977275,  0.46172322,  0.75000097],                         
                          [ 0.59197974,  0.65911348,  0.75000097],                         
                          [ 0.37712092,  0.63661431,  0.75000097],                         
                          [ 0.32023402,  0.46147161,  0.25000057],                         
                          [ 0.48111546,  0.31709408,  0.25000057],                         
                          [ 0.66814176,  0.42548631,  0.25000057],                         
                          [ 0.62284892,  0.63685393,  0.25000057],                         
                          [ 0.40783012,  0.65909408,  0.25000057]]                         
                                                                                           
# Set up configuration                                                                     
bulk_configuration = BulkConfiguration(                                                   
    bravais_lattice=lattice,                                                               
    elements=elements,                                                                     
    fractional_coordinates=fractional_coordinates                                         
    )                                                                                     
                                                                                           
# -------------------------------------------------------------                           
# Calculator                                                                               
# -------------------------------------------------------------                           
#----------------------------------------                                                 
# Basis Set                                                                               
#----------------------------------------                                                 
basis_set = DFTBDirectory("dftb/mio/")                                                     
                                                                                           
#----------------------------------------                                                 
# Pair Potentials                                                                         
#----------------------------------------                                                 
pair_potentials = DFTBDirectory("dftb/mio/")                                               
                                                                                           
numerical_accuracy_parameters = NumericalAccuracyParameters(                               
    k_point_sampling=(1, 1, 10),                                                           
    density_mesh_cutoff=200.0*Hartree,                                                     
    )                                                                                     
                                                                                           
iteration_control_parameters = IterationControlParameters()                               
                                                                                           
poisson_solver = MultigridSolver(                                                         
    boundary_conditions=[[NeumannBoundaryCondition,NeumannBoundaryCondition],             
                         [NeumannBoundaryCondition,NeumannBoundaryCondition],             
                         [PeriodicBoundaryCondition,PeriodicBoundaryCondition]]           
    )                                                                                     
                                                                                           
calculator = SlaterKosterCalculator(                                                       
    basis_set=basis_set,                                                                   
    pair_potentials=pair_potentials,                                                       
    numerical_accuracy_parameters=numerical_accuracy_parameters,                           
    iteration_control_parameters=iteration_control_parameters,                             
    poisson_solver=poisson_solver,                                                         
    )                                                                                     
                                                                                           
bulk_configuration.setCalculator(calculator)                                               
nlprint(bulk_configuration)                                                               
bulk_configuration.update()                                                               
nlsave('F:/P/P-AsP/bs-P-tube-40-dftb-0.2.nc', bulk_configuration)                         
                                                                                           
bulk_configuration = OptimizeGeometry(                                                     
        bulk_configuration,                                                               
        max_forces=0.2*eV/Ang,                                                             
        max_stress=0.2*eV/Ang**3,                                                         
        max_steps=200,                                                                     
        max_step_length=0.2*Ang,                                                           
        trajectory_filename=None,                                                         
        optimizer_method=QuasiNewton(),                                                   
        )                                                                                 
nlsave('F:/P/P-AsP/bs-P-tube-40-dftb-0.2.nc', bulk_configuration)                         
nlprint(bulk_configuration)   
                                                             

[Jess Wellendorff]

I can find nothing wrong with your script as such, so I am guessing the DFTB parameters you use may not be useful for the relaxation. I have attached a script where I relax the forces and the stress (along Z) with DFT-LDA, and that seems to work fine.

[beauyy]

So do you think that the x and y directions should be constrained? Before DFBT, I also used DFT to optimize the structure, similar phenomenon happened.  And I found that the structure suddenly had a great change after max_forces is less than 0.5 eV/Ang. 

[zh]

The vacuum space is included in your simulation unit cell in the a and b (i.e., x and y) directions, so you have to fix them， that is, apply constraint on them. Only the tube length (lattice constant c) and the atomic coordinates are allowed to relax.