Project

General

Profile

CFML_Refcodes.f90

Zaharko, Oksana, 05/02/2012 12:38 PM

Download (360 KB)

 
1
!!----
2
!!---- Copyleft(C) 1999-2011,              Version: 5.0
3
!!---- Juan Rodriguez-Carvajal & Javier Gonzalez-Platas
4
!!----
5
!!---- MODULE: CFML_Keywords_Code_Parser
6
!!----   INFO: Refinable Codes for Parameters
7
!!----
8
!!---- HISTORY
9
!!----    Update: 07/03/2011
10
!!----    Modifications: 12/2011
11
!!----    Modifications: 02/2012
12
!!----
13
!!----
14
!!---- DEPENDENCIES
15
!!----
16
!!----
17
!!---- VARIABLES
18
!!--..    Types
19
!!----    ANGLE_RESTRAINT_TYPE
20
!!----    DISTANCE_RESTRAINT_TYPE
21
!!----    TORSION_RESTRAINT_TYPE
22
!!--..
23
!!----    ANG_REST
24
!!----    CODE_NAM
25
!!----    mCODE_NAM                    !NEW 12/11
26
!!----    DIS_REST
27
!!----    ERR_REFCODES
28
!!----    ERR_REFCODES_MESS
29
!!--++    NCODE                        [Private]
30
!!--++    NKEY                         [Private]
31
!!--++    mNCODE                       [Private] !NEW 12/11
32
!!--++    mNKEY                        [Private] !NEW 12/11
33
!!----    KEY_CODE
34
!!----    KEY_mCODE                    !NEW 12/11
35
!!----    NP_CONS
36
!!----    NP_MAX
37
!!----    NP_REFI
38
!!----    NP_REST_ANG
39
!!----    NP_REST_DIS
40
!!----    NP_REST_TOR
41
!!----    TORSION_REST
42
!!----    V_BCON
43
!!----    V_BOUNDS
44
!!----    V_LIST
45
!!----    V_NAME
46
!!----    V_VEC
47
!!----    V_SHIFT
48
!!----
49
!!---- PUBLIC PROCEDURES
50
!!----    Functions:
51
!!----
52
!!----    Subroutines:
53
!!----       ALLOCATE_RESTPARAM
54
!!----       ALLOCATE_VPARAM
55
!!--++       DELETE_REFCODES             [Private]
56
!!--++       DELETE_REFCODES_FATOM       [Overloaded]
57
!!--++       DELETE_REFCODES_FmATOM      [Overloaded] !NEW 12/11
58
!!--++       DELETE_REFCODES_MOLCRYS     [Overloaded]
59
!!--++       DELETE_REFCODES_MOLEC       [Overloaded]
60
!!--++       DELETE_REFCODES_MAGDOM                   !NEW 02/12
61
!!--++       
62
!!--++       FILL_REFCODES_FATOM         [Overloaded]
63
!!--++       FILL_REFCODES_FmATOM        [Overloaded] !NEW 12/11
64
!!--++       FILL_REFCODES_MOLCRYS       [Overloaded]
65
!!--++       FILL_REFCODES_MOLEC         [Overloaded]
66
!!--++       FILL_REFCODES_MAGDOM                     !NEW 02/12
67
!!--++       GET_ATOMBET_CTR             [Private]
68
!!--++       GET_ATOMPOS_CTR             [Private]
69
!!--++       GET_CONCODES_LINE           [Private]
70
!!--++       GET_CONCODES_LINE_FATOM     [Overloaded]
71
!!--++       GET_CONCODES_LINE_FmATOM    [Overloaded] !NEW 12/11
72
!!--++       GET_CONCODES_LINE_MOLCRYS   [Overloaded]
73
!!--++       GET_CONCODES_LINE_MOLEC     [Overloaded]
74
!!--++       GET_CONCODES_LINE_MAGDOM                 !NEW 02/12
75
!!--++       GET_REFCODES_LINE           [Private]
76
!!--++       GET_REFCODES_LINE_FATOM     [Overloaded]
77
!!--++       GET_REFCODES_LINE_FmATOM    [Overloaded] !NEW 12/11
78
!!--++       GET_REFCODES_LINE_MOLCRYS   [Overloaded]
79
!!--++       GET_REFCODES_LINE_MOLEC     [Overloaded]
80
!!--++       GET_REFCODES_LINE_MAGDOM                 !NEW 02/12
81
!!----       GET_RESTANG_LINE
82
!!----       GET_RESTDIS_LINE
83
!!----       GET_RESTTOR_LINE
84
!!----       INIT_ERR_REFCODES
85
!!----       INIT_REFCODES
86
!!--++       INIT_REFCODES_FATOM         [Overloaded]
87
!!--++       INIT_REFCODES_FmATOM        [Overloaded] !NEW 12/11
88
!!--++       INIT_REFCODES_MOLCRYS       [Overloaded]
89
!!--++       INIT_REFCODES_MOLEC         [Overloaded]
90
!!--++       INIT_REFCODES_MAGDOM                     !NEW 02/12
91
!!----       READ_REFCODES_FILE
92
!!--++       READ_REFCODES_FILE_FATOM    [Overloaded]
93
!!--++       READ_REFCODES_FILE_FmATOM   [Overloaded] !NEW 12/11 !replaced 02/12
94
!!--++       READ_REFCODES_FILE_MagStr                !NEW 02/12
95
!!--++       READ_REFCODES_FILE_MOLCRYS  [Overloaded]
96
!!--++       READ_REFCODES_FILE_MOLEC    [Overloaded]
97
!!--++       SPLIT_OPERATIONS            [Private]
98
!!--++       SPLIT_mOPERATIONS           [Private]    !NEW 12/11
99
!!----       VSTATE_TO_ATOMSPAR
100
!!--++       VSTATE_TO_ATOMSPAR_FATOM    [Overloaded]
101
!!--++       VSTATE_TO_ATOMSPAR_FmATOM   [Overloaded] !NEW 12/11 !modified 02/12 to include MagDom
102
!!--++       VSTATE_TO_ATOMSPAR_MOLCRYS  [Overloaded]
103
!!--++       VSTATE_TO_ATOMSPAR_MOLEC    [Overloaded]
104
!!----       WRITE_INFO_REFCODES
105
!!--++       WRITE_INFO_REFCODES_FATOM   [Overloaded]
106
!!--++       WRITE_INFO_REFCODES_FmATOM  [Overloaded] !NEW 12/11 !replaced 02/12
107
!!--++       WRITE_INFO_REFCODES_Magstr               !NEW 02/12 
108
!!--++       WRITE_INFO_REFCODES_MOLCRYS [Overloaded]
109
!!--++       WRITE_INFO_REFCODES_MOLEC   [Overloaded]
110
!!----       WRITE_INFO_REFPARAMS
111
!!----       WRITE_RESTRAINTS_OBSCALC
112
!!----
113
!!
114
 Module CFML_Keywords_Code_Parser
115
    !---- Modules ----!
116
    Use CFML_GlobalDeps,                only: cp
117
    Use CFML_Math_General,              only: Sort
118
    Use CFML_String_Utilities,          only: Cutst, U_Case, L_Case, Getword, GetNum
119
    Use CFML_Crystallographic_Symmetry, only: Space_Group_Type, Get_Stabilizer, Symmetry_Symbol,   &
120
                                              Sym_B_Relations, Read_SymTrans_Code, Get_SymSymb,Read_Xsym
121
    Use CFML_Atom_TypeDef,              only: Atom_list_Type, mAtom_list_Type
122
    Use CFML_Molecular_Crystals,        only: Molecule_Type, Molecular_Crystal_Type
123
    Use CFML_IO_Formats,                only: File_List_Type
124
    Use CFML_Magnetic_Symmetry
125
    Use CFML_Math_3D,                   only: Get_Spheric_Coord
126

    
127
    !---- Variables ----!
128
    implicit none
129

    
130
    private
131

    
132
    !---- List of public functions ----!
133

    
134
    !---- List of public subroutines ----!
135
    public :: Allocate_VParam, Init_RefCodes, Read_RefCodes_File, VState_to_AtomsPar,  &
136
              Write_Info_RefCodes, Get_RestAng_Line, Get_RestDis_Line, Get_RestTor_Line, &
137
              Allocate_RestParam, Write_Restraints_ObsCalc, Init_Err_RefCodes, &
138
              Write_Info_RefParams
139

    
140
    !---- List of private functions ----!
141

    
142
    !---- List of private subroutines ----!
143
    private :: Delete_RefCodes,   &
144
               Delete_RefCodes_FAtom, Delete_RefCodes_FmAtom, Delete_RefCodes_Molcrys,          &
145
               Delete_RefCodes_Molec, Delete_RefCodes_Magdom, &
146
               Fill_RefCodes,     &
147
               Fill_RefCodes_FAtom, Fill_RefCodes_FmAtom, Fill_RefCodes_Molcrys,                &
148
               Fill_RefCodes_Molec, Fill_RefCodes_Magdom,  &
149
               Get_AtomBet_Ctr, Get_Atompos_Ctr,                                                &
150
               Get_ConCodes_Line, &
151
               Get_ConCodes_Line_FAtom, Get_ConCodes_Line_FmAtom, Get_ConCodes_Line_Molcrys,    &
152
               Get_ConCodes_Line_Molec, Get_ConCodes_Line_Magdom, &
153
               Get_RefCodes_Line, &
154
               Get_RefCodes_Line_FAtom, Get_RefCodes_Line_FmAtom, Get_RefCodes_Line_Molcrys,    &
155
               Get_RefCodes_Line_Molec, Get_RefCodes_Line_Magdom, &
156
               Init_RefCodes_FAtom, Init_RefCodes_FmAtom, Init_RefCodes_Molcrys,                &
157
               Init_RefCodes_Molec, Init_RefCodes_Magdom, &
158
               Read_RefCodes_File_FAtom, Read_RefCodes_File_MagStr, Read_RefCodes_File_Molcrys, &
159
               Read_RefCodes_File_Molec, &
160
               Split_Operations, Split_mOperations, &
161
               VState_to_AtomsPar_FAtom, VState_to_AtomsPar_FmAtom, VState_to_AtomsPar_Molcrys, &
162
               VState_to_AtomsPar_Molec, &
163
               Write_Info_RefCodes_FAtom,Write_Info_RefCodes_MagStr,Write_Info_RefCodes_Molcrys,& 
164
               Write_Info_RefCodes_Molec
165

    
166
    !---- Definitions ----!
167

    
168
    !!----
169
    !!---- TYPE :: ANGLE_RESTRAINT_TYPE
170
    !!--..
171
    !!---- Type, public :: Angle_Restraint_Type
172
    !!----    real(kind=cp)                 :: AObs
173
    !!----    real(kind=cp)                 :: ACalc
174
    !!----    real(kind=cp)                 :: Sigma
175
    !!----    integer,dimension(3)          :: P
176
    !!----    character(len=8),dimension(2) :: STCode
177
    !!---- End Type Angle_Restraint_Type
178
    !!----
179
    !!---- Update: April - 2005
180
    !!
181
    Type, public :: Angle_Restraint_Type
182
      real(kind=cp)                 :: AObs
183
      real(kind=cp)                 :: ACalc
184
      real(kind=cp)                 :: Sigma
185
      integer,dimension(3)          :: P
186
      character(len=8),dimension(2) :: STCode
187
    End Type Angle_Restraint_Type
188

    
189
    !!----
190
    !!---- TYPE :: DISTANCE_RESTRAINT_TYPE
191
    !!--..
192
    !!---- Type, public :: Distance_Restraint_Type
193
    !!----    real(kind=cp)        :: DObs
194
    !!----    real(kind=cp)        :: DCalc
195
    !!----    real(kind=cp)        :: Sigma
196
    !!----    integer,dimension(2) :: P
197
    !!----    character(len=8)     :: STCode    ! _N.ABC
198
    !!---- End Type Distance_Restraint_Type
199
    !!----
200
    !!---- Update: April - 2005
201
    !!
202
    Type, public :: Distance_Restraint_Type
203
      real(kind=cp)        :: DObs
204
      real(kind=cp)        :: DCalc
205
      real(kind=cp)        :: Sigma
206
      integer,dimension(2) :: P
207
      character(len=8)     :: STCode
208
    End Type Distance_Restraint_Type
209

    
210
    !!----
211
    !!---- TYPE :: TORSION_RESTRAINT_TYPE
212
    !!--..
213
    !!---- Type, public :: Torsion_Restraint_Type
214
    !!----    real(kind=cp)                 :: TObs
215
    !!----    real(kind=cp)                 :: TCalc
216
    !!----    real(kind=cp)                 :: Sigma
217
    !!----    integer,dimension(4)          :: P
218
    !!----    character(len=8),dimension(3) :: STCode
219
    !!---- End Type Torsion_Restraint_Type
220
    !!----
221
    !!---- Update: April - 2005
222
    !!
223
    Type, public :: Torsion_Restraint_Type
224
      real(kind=cp)                 :: TObs
225
      real(kind=cp)                 :: TCalc
226
      real(kind=cp)                 :: Sigma
227
      integer,dimension(4)          :: P
228
      character(len=8),dimension(3) :: STCode
229
    End Type Torsion_Restraint_Type
230

    
231
    !!----
232
    !!---- ANG_REST
233
    !!----    type(Angle_Restraint_Type), public, dimension(:), allocatable :: Ang_rest
234
    !!----
235
    !!---- Relations for Angle Restraints
236
    !!----
237
    !!---- Update: March - 2005
238
    !!
239
    type(Angle_Restraint_Type), public, dimension(:), allocatable :: Ang_Rest
240

    
241
    !!--++
242
    !!--++ NCODE
243
    !!--++    integer, private :: NCode
244
    !!--++
245
    !!--++    Number of Code variables
246
    !!--++
247
    !!--++ Update: March - 2005
248
    !!
249
    integer, private, parameter :: NCode=21
250

    
251
    !!----
252
    !!---- CODE_NAM
253
    !!----    character(len=*), dimension(NCode), public, parameter :: Code_Nam
254
    !!----
255
    !!----    Variable for treatement Codes
256
    !!--..    21 debe ser igual a NCode
257
    !!----
258
    !!---- Update: March - 2005
259
    !!
260
    character(len=*), dimension(NCode), public, parameter :: Code_Nam=(/ "X_    ","Y_    ","Z_    ",       &
261
                                                                         "B_    ","Occ_  ","B11_  ",       &
262
                                                                         "B22_  ","B33_  ","B12_  ",       &
263
                                                                         "B13_  ","B23_  ","Bns_  ",       &
264
                                                                         "Xc_   ","Yc_   ","Zc_   ",       &
265
                                                                         "Theta_","Phi_  ","Chi_  ",       &
266
                                                                         "Th_L_ ","Th_T_ ","Th_S_ "/)
267
    integer, private, parameter :: mNCode=25
268

    
269
    !!----
270
    !!---- mCODE_NAM
271
    !!----    character(len=*), dimension(mNCode), public, parameter :: mCode_Nam
272
    !!----
273
    !!----    Variable for treatement Codes
274
    !!----    mag clone of CODE_NAM
275
    !!---- Created: December - 2011
276
    !!
277
    character(len=*),dimension(mNCode), public, parameter :: mCode_Nam=(/"Rx_   ","Ry_   ","Rz_   ",&
278
                                                                         "Ix_   ","Iy_   ","Iz_   ",&
279
                                                                         "Rm_   ","Rphi_ ","Rth_  ",&
280
                                                                         "Im_   ","Iphi_ ","Ith_  ",&
281
                                                                         "MagPh_",                  &
282
                                                                         "C1_   ","C2_   ","C3_   ",&
283
                                                                         "C4_   ","C5_   ","C6_   ",&
284
                                                                         "C7_   ","C8_   ","C9_   ",&
285
                                                                         "C10_  ","C11_  ","C12_  "/)
286

    
287
    !!----
288
    !!---- DIS_REST
289
    !!----    type(Distance_Restraint_Type), public, dimension(:), allocatable :: Dis_Rest
290
    !!----
291
    !!---- Relations for Angle Restraints
292
    !!----
293
    !!---- Update: March - 2005
294
    !!
295
    type(Distance_Restraint_Type), public, dimension(:), allocatable :: Dis_Rest
296

    
297
    !!----
298
    !!---- ERR_REFCODES
299
    !!----    logical, public :: Err_RefCodes
300
    !!----
301
    !!----    Error variable
302
    !!----
303
    !!---- Update: March - 2005
304
    !!
305
    logical, public :: Err_RefCodes = .false.
306

    
307
    !!----
308
    !!---- ERR_REFCODES_MESS
309
    !!----    character(len=150), public :: ERR_RefCodes_Mess
310
    !!----
311
    !!----    Error variable messages
312
    !!----
313
    !!---- Update: March - 2005
314
    !!
315
    character(len=150), public :: ERR_RefCodes_Mess = " "
316

    
317
    !!--++
318
    !!--++ NKEY
319
    !!--++    integer, public :: NKey
320
    !!--++
321
    !!--++    Number of Keys variables
322
    !!--++
323
    !!--++ Update: March - 2005
324
    !!
325
    integer, private, parameter :: NKey=8
326

    
327
    !!----
328
    !!---- KEY_CODE
329
    !!----    character(len=*), dimension(NKey), public, parameter :: key_Code
330
    !!----
331
    !!----     Key codes defined in the module
332
    !!----
333
    !!---- Update: March - 2005
334
    !!
335
    character(len=*), dimension(Nkey), public, parameter :: Key_Code=(/ "XYZ ","OCC ","BIS ","BAN ", &
336
                                                                        "ALL ","CEN ","ORI ","THE "/)
337
    integer, private, parameter :: mNKey=4
338

    
339
    !!----
340
    !!---- KEY_mCODE
341
    !!----    character(len=*), dimension(mNKey), public, parameter :: key_mCode
342
    !!----
343
    !!----    Key codes defined in the module
344
    !!----    mag clone of KEY_CODE
345
    !!---- Created: December - 2011
346
    !!
347
    character(len=*), dimension(mNkey), public, parameter :: Key_mCode=(/"Rxyz","Ixyz","Mxyz","Magd"/)
348

    
349
    !!----
350
    !!---- NP_CONS
351
    !!----    integer, public :: NP_Cons
352
    !!----
353
    !!----    Number of Constraints relations
354
    !!----
355
    !!---- Update: March - 2005
356
    !!
357
    integer, public :: NP_Cons
358

    
359
    !!----
360
    !!---- NP_MAX
361
    !!----    integer, public :: NP_Max
362
    !!----
363
    !!----    Number of Maximum Parameters to Refine
364
    !!----
365
    !!---- Update: March - 2005
366
    !!
367
    integer, public :: NP_Max
368

    
369
    !!----
370
    !!---- NP_REFI
371
    !!----    integer, public :: NP_Refi
372
    !!----
373
    !!----    Number of Refinable Parameters
374
    !!----
375
    !!---- Update: March - 2005
376
    !!
377
    integer, public :: NP_Refi
378

    
379
    !!----
380
    !!---- NP_REST_ANG
381
    !!----    integer, public :: NP_Rest_Ang
382
    !!----
383
    !!----    Number of Angle Restraints relations
384
    !!----
385
    !!---- Update: March - 2005
386
    !!
387
    integer, public :: NP_Rest_Ang
388

    
389
    !!----
390
    !!---- NP_REST_DIS
391
    !!----    integer, public :: NP_Rest_Dis
392
    !!----
393
    !!----    Number of Distance Restraints relations
394
    !!----
395
    !!---- Update: March - 2005
396
    !!
397
    integer, public :: NP_Rest_Dis
398

    
399
    !!----
400
    !!---- NP_REST_TOR
401
    !!----    integer, public :: NP_Rest_Tor
402
    !!----
403
    !!----    Number of Torsion Restraints relations
404
    !!----
405
    !!---- Update: March - 2005
406
    !!
407
    integer, public :: NP_Rest_Tor
408

    
409
    !!----
410
    !!---- TOR_REST
411
    !!----    type(Torsion_Restraint_Type), public, dimension(:), allocatable :: Tor_Rest
412
    !!----
413
    !!---- Relations for Torsion Angle Restraints
414
    !!----
415
    !!---- Update: March - 2005
416
    !!
417
    type(Torsion_Restraint_Type), public, dimension(:), allocatable :: Tor_Rest
418

    
419
    !!----
420
    !!---- V_BCON
421
    !!----    integer, public, dimension(:), allocatable :: V_BCon
422
    !!----
423
    !!----    Vector of Boundary Conditions
424
    !!----
425
    !!---- Update: March - 2005
426
    !!
427
    integer, public, dimension(:), allocatable :: V_BCon
428

    
429
    !!----
430
    !!---- V_BOUNDS
431
    !!----    real(kind=cp), public, dimension(:,:), allocatable :: V_Bounds
432
    !!----
433
    !!----    Vector of Lower, Upper limits and Step for Parameters
434
    !!----
435
    !!---- Update: March - 2005
436
    !!
437
    real(kind=cp), public, dimension(:,:),  allocatable :: V_Bounds
438

    
439
    !!----
440
    !!---- V_LIST
441
    !!----    integer, public, dimension(:), allocatable :: V_List
442
    !!----
443
    !!----    Vector of Index point the atom order
444
    !!----
445
    !!---- Update: March - 2005
446
    !!
447
    integer, public, dimension(:),  allocatable :: V_List
448

    
449
    !!----
450
    !!---- V_NAME
451
    !!----    character(len=20), public, dimension(:), allocatable :: V_Name
452
    !!----
453
    !!----    Vector of  Name of Refinable Parameters
454
    !!----
455
    !!---- Update: March - 2005
456
    !!
457
    character(len=20), public, dimension(:), allocatable :: V_Name
458

    
459
    !!----
460
    !!---- V_VEC
461
    !!----    real(kind=cp), public, dimension(:), allocatable :: V_Vec
462
    !!----
463
    !!----    Vector of  Parameters
464
    !!----
465
    !!---- Update: March - 2005
466
    !!
467
    real(kind=cp), public, dimension(:),    allocatable :: V_Vec
468

    
469
    !!----
470
    !!---- V_SHIFT
471
    !!----    real(kind=cp), public, dimension(:), allocatable :: V_Shift
472
    !!----
473
    !!----    Vector of holding the shift of parameters
474
    !!----
475
    !!---- Update: March - 2005
476
    !!
477
    real(kind=cp), public, dimension(:),    allocatable :: V_Shift
478

    
479
    !---- Interfaces - Overloaded ----!
480
    Interface Delete_RefCodes
481
       Module Procedure Delete_RefCodes_FAtom
482
       Module Procedure Delete_RefCodes_FmAtom
483
       Module Procedure Delete_RefCodes_Molcrys
484
       Module Procedure Delete_RefCodes_Molec
485
       Module Procedure Delete_RefCodes_Magdom
486
    End Interface
487

    
488
    Interface Fill_RefCodes
489
       Module Procedure Fill_RefCodes_FAtom
490
       Module Procedure Fill_RefCodes_FmAtom
491
       Module Procedure Fill_RefCodes_Molcrys
492
       Module Procedure Fill_RefCodes_Molec
493
       Module Procedure Fill_RefCodes_Magdom
494
    End Interface
495

    
496
    Interface Get_ConCodes_Line
497
       Module Procedure Get_ConCodes_Line_FAtom
498
       Module Procedure Get_ConCodes_Line_FmAtom
499
       Module Procedure Get_ConCodes_Line_Molcrys
500
       Module Procedure Get_ConCodes_Line_Molec
501
       Module Procedure Get_ConCodes_Line_Magdom
502
    End Interface
503

    
504
    Interface Get_RefCodes_Line
505
       Module Procedure Get_RefCodes_Line_FAtom
506
       Module Procedure Get_RefCodes_Line_FmAtom
507
       Module Procedure Get_RefCodes_Line_Molcrys
508
       Module Procedure Get_RefCodes_Line_Molec
509
       Module Procedure Get_RefCodes_Line_Magdom
510
    End Interface
511

    
512
    Interface Init_RefCodes
513
       Module Procedure Init_RefCodes_FAtom
514
       Module Procedure Init_RefCodes_FmAtom
515
       Module Procedure Init_RefCodes_Molcrys
516
       Module Procedure Init_RefCodes_Molec
517
       Module Procedure Init_RefCodes_Magdom
518
    End Interface
519

    
520
    Interface Read_RefCodes_File
521
       Module Procedure Read_RefCodes_File_FAtom
522
       Module Procedure Read_RefCodes_File_MagStr
523
       Module Procedure Read_RefCodes_File_Molcrys
524
       Module Procedure Read_RefCodes_File_Molec
525
    End Interface
526

    
527
    Interface VState_to_AtomsPar
528
       Module Procedure VState_to_AtomsPar_FAtom
529
       Module Procedure VState_to_AtomsPar_FmAtom
530
       Module Procedure VState_to_AtomsPar_Molcrys
531
       Module Procedure VState_to_AtomsPar_Molec
532
    End Interface
533

    
534
    Interface Write_Info_RefCodes
535
       Module Procedure Write_Info_RefCodes_FAtom
536
       Module Procedure Write_Info_RefCodes_MagStr
537
       Module Procedure Write_Info_RefCodes_Molcrys
538
       Module Procedure Write_Info_RefCodes_Molec
539
    End Interface
540

    
541
 Contains
542

    
543
    !!----
544
    !!---- Subroutine Allocate_RestParam(file_dat)
545
    !!----    Type(file_list_type),     intent( in)    :: file_dat
546
    !!----
547
    !!----    Allocate vectors Ang_Rest, Dist_Rest, Tor_Rest
548
    !!----
549
    !!---- Update: March - 2005
550
    !!
551
    Subroutine Allocate_RestParam(file_dat)
552
       !---- Arguments ----!
553
       Type(file_list_type),     intent( in)    :: file_dat
554

    
555
       !---- Local variables ----!
556
       character(len=132)              :: line
557
       character(len=15),dimension(40) :: car
558
       integer                         :: i,nc,nr
559

    
560
       if (allocated(Ang_Rest)) deallocate(Ang_Rest)
561
       if (allocated(Dis_Rest)) deallocate(Dis_Rest)
562
       if (allocated(Tor_Rest)) deallocate(Tor_Rest)
563

    
564
       NP_Rest_Ang=0
565
       NP_Rest_Dis=0
566
       NP_Rest_Tor=0
567

    
568
       !---- Dimension for AFIX ----!
569
       nr=0
570
       do i=1,file_dat%nlines
571
          line=adjustl(file_dat%line(i))
572
          if (u_case(line(1:4)) /= "AFIX") cycle
573
          call cutst(line)
574
          call getword(line,car,nc)
575
          nr=nr+nc/3
576
       end do
577
       if (nr >0) then
578
          allocate(Ang_Rest(nr))
579
          ang_rest%aobs =0.0
580
          ang_rest%acalc=0.0
581
          ang_rest%sigma=0.0
582
          ang_rest%p(1) = 0
583
          ang_rest%p(2) = 0
584
          ang_rest%STCode(1)=" "
585
          ang_rest%STCode(2)=" "
586
       end if
587

    
588
       !---- Dimension for DFIX ----!
589
       nr=0
590
       do i=1,file_dat%nlines
591
          line=adjustl(file_dat%line(i))
592
          if (u_case(line(1:4)) /= "DFIX") cycle
593
          call cutst(line)
594
          call getword(line,car,nc)
595
          nr=nr+nc/2
596
          if (modulo(nc,2) == 0) nr=nr-1
597
       end do
598
       if (nr >0) then
599
          allocate(Dis_Rest(nr))
600
          dis_rest%dobs =0.0
601
          dis_rest%dcalc=0.0
602
          dis_rest%sigma=0.0
603
          dis_rest%p(1) = 0
604
          dis_rest%p(2) = 0
605
          dis_rest%STCode=" "
606
       end if
607

    
608
       !---- Dimension for TFIX ----!
609
       nr=0
610
       do i=1,file_dat%nlines
611
          line=adjustl(file_dat%line(i))
612
          if (u_case(line(1:4)) /= "TFIX") cycle
613
          call cutst(line)
614
          call getword(line,car,nc)
615
          nr=nr+nc/4
616
       end do
617
       if (nr >0) then
618
          allocate(Tor_Rest(nr))
619
          tor_rest%tobs=0.0
620
          tor_rest%tcalc=0.0
621
          tor_rest%sigma=0.0
622
          tor_rest%p(1) = 0
623
          tor_rest%p(2) = 0
624
          tor_rest%STCode(1)=" "
625
          tor_rest%STCode(2)=" "
626
          tor_rest%STCode(3)=" "
627
       end if
628

    
629
       return
630
    End Subroutine Allocate_RestParam
631

    
632
    !!---- Subroutine Allocate_VParam(N)
633
    !!----    integer, intent(in) :: N
634
    !!----
635
    !!----    Allocate vectors V_Vec, V_Bounds, V_Name, V_Bcon, V_Shift, V_list
636
    !!----    If N is equal zero it deallocates the vectors
637
    !!----
638
    !!---- Update: March - 2005
639
    !!
640
    Subroutine Allocate_VParam(N)
641
       !---- Arguments ----!
642
       integer, intent(in) :: N
643

    
644
       if (allocated(V_Vec))    deallocate(V_Vec)
645
       if (allocated(V_Name))   deallocate(V_Name)
646
       if (allocated(V_Bounds)) deallocate(V_Bounds)
647
       if (allocated(V_BCon))   deallocate(V_BCon)
648
       if (allocated(V_Shift))  deallocate(V_Shift)
649
       if (allocated(V_List))   deallocate(V_List)
650

    
651
       if (N > 0) then
652
          allocate(V_Vec(n))
653
          V_Vec=0.0
654
          allocate(V_Name(n))
655
          V_Name=" "
656
          allocate(V_Bounds(3,n))
657
          V_Bounds=0.0
658
          allocate(V_BCon(n))
659
          V_BCon=0
660
          allocate(V_Shift(n))
661
          V_Shift=0.0
662
          allocate(V_List(n))
663
          V_list=0
664
          np_max=n
665
       else
666
          np_max=0
667
       end if
668

    
669
       return
670
    End Subroutine Allocate_VParam
671

    
672
    !!--++
673
    !!--++ Subroutine Delete_RefCodes(N, FAtom/FmAtom/MolCrys/Molec/Mag_Dom)
674
    !!--++    integer,                      intent(in)     :: N
675
    !!--++    type(Atom_List_Type),         intent(in out) :: FAtom
676
    !!--++    or
677
    !!--++    type(mAtom_List_Type),        intent(in out) :: FmAtom
678
    !!--++    or
679
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
680
    !!--++    or
681
    !!--++    type(molecule_type),          intent(in out) :: Molec
682
    !!--++    or
683
    !!--++    type(Magnetic_Domain_type),intent(in out) :: Mag_Dom
684
    !!--++
685
    !!--++    (Private)
686
    !!--++    Delete the number of Refinable Parameter (N) on the list
687
    !!--++
688
    !!--++ Update: March - 2005
689
    !!
690

    
691
    !!--++
692
    !!--++ Subroutine Delete_RefCodes_FAtom(N, FAtom)
693
    !!--++    integer,              intent(in)     :: N
694
    !!--++    type(Atom_List_Type), intent(in out) :: FAtom
695
    !!--++
696
    !!--++ Overloaded
697
    !!--++ Delete the number of Refinable Parameter (N) on the list
698
    !!--++
699
    !!--++ Update: March - 2005
700
    !!
701
    Subroutine Delete_RefCodes_FAtom(N, FAtom)
702
       !---- Arguments ----!
703
       integer,              intent(in)     :: N
704
       type(Atom_List_Type), intent(in out) :: FAtom
705

    
706
       !---- Local Variables ----!
707
       logical :: deleted
708
       integer :: i,j
709

    
710
       deleted=.false.
711

    
712
       !---- Eliminate N Parameter ----!
713
       do i=1,FAtom%natoms
714
          do j=1,3
715
             if (FAtom%atom(i)%lx(j) == N) then
716
                 FAtom%atom(i)%lx(j)=0
717
                 FAtom%atom(i)%mx(j)=0.0
718
                 deleted=.true.
719
             end if
720
          end do
721

    
722
          if (FAtom%atom(i)%lbiso == N) then
723
              FAtom%atom(i)%lbiso=0
724
              FAtom%atom(i)%mbiso=0.0
725
              deleted=.true.
726
          end if
727

    
728
          if (FAtom%atom(i)%locc == N) then
729
              FAtom%atom(i)%locc=0
730
              FAtom%atom(i)%mocc=0.0
731
              deleted=.true.
732
          end if
733

    
734
          do j=1,6
735
             if (FAtom%atom(i)%lu(j) == N) then
736
                 FAtom%atom(i)%lu(j)=0
737
                 FAtom%atom(i)%mu(j)=0.0
738
                 deleted=.true.
739
             end if
740
          end do
741
       end do
742

    
743
       !---- Updating Variables ----!
744
       do i=1,FAtom%natoms
745
          do j=1,3
746
             if (FAtom%atom(i)%lx(j) > N) then
747
                FAtom%atom(i)%lx(j)=FAtom%atom(i)%lx(j)-1
748
             end if
749
          end do
750

    
751
          if (FAtom%atom(i)%lbiso > N) then
752
             FAtom%atom(i)%lbiso=FAtom%atom(i)%lbiso-1
753
          end if
754

    
755
          if (FAtom%atom(i)%locc > N) then
756
             FAtom%atom(i)%locc=FAtom%atom(i)%locc-1
757
          end if
758

    
759
          do j=1,6
760
             if (FAtom%atom(i)%lu(j) > N) then
761
                FAtom%atom(i)%lu(j)=FAtom%atom(i)%lu(j)-1
762
             end if
763
          end do
764
       end do
765

    
766
       !---- Updating V_Vectors ----!
767
       if (deleted) then
768
          do i=N+1,Np_refi
769
             V_Vec(i-1)=V_Vec(i)
770
             V_Name(i-1)=V_Name(i)
771
             V_Bounds(:,i-1)=V_Bounds(:,i)
772
             V_BCon(i-1)=V_BCon(i)
773
             V_List(i-1)=V_List(i)
774
          end do
775
          V_Vec(np_refi)=0.0
776
          V_Name(np_refi)=" "
777
          V_Bounds(:,np_refi)=0.0
778
          V_BCon(np_refi)=0
779
          V_List(np_refi)=0
780

    
781
          np_refi=np_refi-1
782
       end if
783

    
784
       return
785
    End Subroutine Delete_RefCodes_FAtom
786

    
787
    !!--++
788
    !!--++ Subroutine Delete_RefCodes_FmAtom(N, FmAtom)
789
    !!--++    integer,              intent(in)     :: N
790
    !!--++    type(mAtom_List_Type),intent(in out) :: FmAtom
791
    !!--++
792
    !!--++ Delete the number of Refinable Parameters (N) on the list
793
    !!--++ magnetic clone of Delete_RefCodes_FAtom
794
    !!--++ Created: December - 2011
795
    !!
796
    Subroutine Delete_RefCodes_FmAtom(N, FmAtom)
797
       !---- Arguments ----!
798
       integer,              intent(in)     :: N
799
       type(mAtom_List_Type), intent(in out) :: FmAtom
800

    
801
       !---- Local Variables ----!
802
       logical :: deleted
803
       integer :: i,j,k,ik
804

    
805
       deleted=.false.
806

    
807
       !---- Eliminate N Parameter ----!
808
       do i=1,FmAtom%natoms
809
       ik=FmAtom%atom(i)%nvk
810
          do j=1,3
811
           do k=1,ik
812
             if (FmAtom%atom(i)%lSkR(j,k) == N) then
813
                 FmAtom%atom(i)%mSkR(j,k)=0.0
814
                 FmAtom%atom(i)%lskr(j,k)=0
815
                 deleted=.true.
816
             end if
817
           end do
818
          end do
819

    
820
          do j=1,3
821
           do k=1,ik
822
             if (FmAtom%atom(i)%lSkI(j,k) == N) then
823
                 FmAtom%atom(i)%mSkI(j,k)=0.0
824
                 FmAtom%atom(i)%lski(j,k)=0
825
                 deleted=.true.
826
             end if
827
           end do
828
          end do
829

    
830
          do k=1,ik
831
             if (FmAtom%atom(i)%lmphas(k) == N) then
832
                 FmAtom%atom(i)%mmphas(k)=0.0
833
                 FmAtom%atom(i)%lmphas(k)=0
834
                 deleted=.true.
835
             end if
836
          end do
837

    
838
          do j=1,12
839
           do k=1,ik
840
             if (FmAtom%atom(i)%lbas(j,k) == N) then
841
                 FmAtom%atom(i)%mbas(j,k)=0.0
842
                 FmAtom%atom(i)%lbas(j,k)=0
843
                 deleted=.true.
844
             end if
845
           end do
846
          end do
847
       end do
848

    
849
       !---- Updating Variables ----!
850
       do i=1,FmAtom%natoms
851
          do j=1,3
852
           do k=1,ik
853
             if (FmAtom%atom(i)%lSkR(j,k) > N) then
854
                 FmAtom%atom(i)%lSkR(j,k)=FmAtom%atom(i)%lSkR(j,k)-1
855
             end if
856
           end do
857
          end do
858
          
859
          do j=1,3
860
           do k=1,ik
861
             if (FmAtom%atom(i)%lSkI(j,k) > N) then
862
                 FmAtom%atom(i)%lSkI(j,k)=FmAtom%atom(i)%lSkI(j,k)-1
863
             end if
864
           end do
865
          end do
866

    
867
          do k=1,ik
868
             if (FmAtom%atom(i)%lmphas(k) > N) then
869
                 FmAtom%atom(i)%lmphas(k)=FmAtom%atom(i)%lmphas(k)-1
870
             end if
871
          end do
872

    
873
          do j=1,12
874
           do k=1,ik
875
             if (FmAtom%atom(i)%lbas(j,k) > N) then
876
                 FmAtom%atom(i)%lbas(j,k)=FmAtom%atom(i)%lbas(j,k)-1
877
             end if
878
           end do
879
          end do
880
       end do
881

    
882
       !---- Updating V_Vectors ----!
883
       if (deleted) then
884
          do i=N+1,Np_refi
885
             V_Vec(i-1)=V_Vec(i)
886
             V_Name(i-1)=V_Name(i)
887
             V_Bounds(:,i-1)=V_Bounds(:,i)
888
             V_BCon(i-1)=V_BCon(i)
889
             V_List(i-1)=V_List(i)
890
          end do
891
          V_Vec(np_refi)=0.0
892
          V_Name(np_refi)=" "
893
          V_Bounds(:,np_refi)=0.0
894
          V_BCon(np_refi)=0
895
          V_List(np_refi)=0
896

    
897
          np_refi=np_refi-1
898
       end if
899

    
900
       return
901
    End Subroutine Delete_RefCodes_FmAtom
902

    
903
    !!--++
904
    !!--++ Subroutine Delete_RefCodes_MolCrys(N,MolCrys)
905
    !!--++    integer,                      intent(in)     :: N
906
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
907
    !!--++
908
    !!--++ Overloaded
909
    !!--++ Delete the number of Refinable Parameter (N) on the list
910
    !!--++
911
    !!--++ Update: March - 2005
912
    !!
913
    Subroutine Delete_RefCodes_MolCrys(N,Molcrys)
914
       !---- Arguments ----!
915
       integer,                      intent(in)     :: N
916
       type(molecular_Crystal_type), intent(in out) :: MolCrys
917

    
918
       !---- Local Variables ----!
919
       logical :: deleted
920
       integer :: i,j,k
921

    
922
       deleted=.false.
923

    
924
       if (MolCrys%N_Free > 0 ) then
925
          do i=1,MolCrys%N_Free
926
             do j=1,3
927
                if (MolCrys%Atm(i)%lx(j) == N) then
928
                    MolCrys%Atm(i)%lx(j)=0
929
                    MolCrys%Atm(i)%mx(j)=0.0
930
                    deleted=.true.
931
                end if
932
             end do
933

    
934
             if (MolCrys%Atm(i)%lbiso == N) then
935
                 MolCrys%Atm(i)%lbiso=0
936
                 MolCrys%Atm(i)%mbiso=0.0
937
                 deleted=.true.
938
             end if
939

    
940
             if (MolCrys%Atm(i)%locc == N) then
941
                 MolCrys%Atm(i)%locc=0
942
                 MolCrys%Atm(i)%mocc=0.0
943
                 deleted=.true.
944
             end if
945

    
946
             do j=1,6
947
                if (MolCrys%Atm(i)%lu(j) == N) then
948
                    MolCrys%Atm(i)%lu(j)=0
949
                    MolCrys%Atm(i)%mu(j)=0.0
950
                    deleted=.true.
951
                end if
952
             end do
953
          end do
954

    
955
          do i=1,MolCrys%N_Free
956
             do j=1,3
957
                if (MolCrys%Atm(i)%lx(j) > N) then
958
                    MolCrys%Atm(i)%lx(j)=MolCrys%Atm(i)%lx(j)-1
959
                end if
960
             end do
961

    
962
             if (MolCrys%Atm(i)%lbiso > N) then
963
                MolCrys%Atm(i)%lbiso=MolCrys%Atm(i)%lbiso-1
964
             end if
965

    
966
             if (MolCrys%Atm(i)%locc > N) then
967
                MolCrys%Atm(i)%locc=MolCrys%Atm(i)%locc-1
968
             end if
969

    
970
             do j=1,6
971
                if (MolCrys%Atm(i)%lu(j) > N) then
972
                   MolCrys%Atm(i)%lu(j)=MolCrys%Atm(i)%lu(j)-1
973
                end if
974
             end do
975
          end do
976
       end if
977

    
978
       if (MolCrys%N_Mol > 0 ) then
979

    
980
          do k=1,MolCrys%N_Mol
981
             do j=1,3
982
                if (Molcrys%Mol(k)%lxcentre(j) == N) then
983
                   Molcrys%Mol(k)%lxcentre(j)=0
984
                   Molcrys%Mol(k)%mxcentre(j)=0.0
985
                   deleted=.true.
986
                end if
987
             end do
988

    
989
             do j=1,3
990
                if (Molcrys%Mol(k)%lOrient(j) == N) then
991
                   Molcrys%Mol(k)%lOrient(j)=0
992
                   Molcrys%Mol(k)%mOrient(j)=0.0
993
                   deleted=.true.
994
                end if
995
             end do
996

    
997
             do j=1,6
998
                if (Molcrys%Mol(k)%lT_TLS(j) == N) then
999
                   Molcrys%Mol(k)%lT_TLS(j)=0
1000
                   Molcrys%Mol(k)%mT_TLS(j)=0.0
1001
                   deleted=.true.
1002
                end if
1003
             end do
1004

    
1005
             do j=1,6
1006
                if (Molcrys%Mol(k)%lL_TLS(j) == N) then
1007
                   Molcrys%Mol(k)%lL_TLS(j)=0
1008
                   Molcrys%Mol(k)%mL_TLS(j)=0.0
1009
                   deleted=.true.
1010
                end if
1011
             end do
1012

    
1013
             do i=1,3
1014
                do j=1,3
1015
                   if (Molcrys%Mol(k)%lS_TLS(i,j) == N) then
1016
                      Molcrys%Mol(k)%lS_TLS(i,j)=0
1017
                      Molcrys%Mol(k)%mS_TLS(i,j)=0.0
1018
                      deleted=.true.
1019
                   end if
1020
                end do
1021
             end do
1022

    
1023
             !---- Updating ----!
1024
             do j=1,3
1025
                if (Molcrys%Mol(k)%lxcentre(j) > N) then
1026
                   Molcrys%Mol(k)%lxcentre(j)=Molcrys%Mol(k)%lxcentre(j)-1
1027
                end if
1028
             end do
1029

    
1030
             do j=1,3
1031
                if (Molcrys%Mol(k)%lOrient(j) > N) then
1032
                   Molcrys%Mol(k)%lOrient(j)=Molcrys%Mol(k)%lOrient(j)-1
1033
                end if
1034
             end do
1035

    
1036
             do j=1,6
1037
                if (Molcrys%Mol(k)%lT_TLS(j) > N) then
1038
                   Molcrys%Mol(k)%lT_TLS(j)=Molcrys%Mol(k)%lT_TLS(j)-1
1039
                end if
1040
             end do
1041

    
1042
             do j=1,6
1043
                if (Molcrys%Mol(k)%lL_TLS(j) > N) then
1044
                   Molcrys%Mol(k)%lL_TLS(j)=Molcrys%Mol(k)%lL_TLS(j)-1
1045
                end if
1046
             end do
1047

    
1048
             do i=1,3
1049
                do j=1,3
1050
                   if (Molcrys%Mol(k)%lS_TLS(i,j) > N) then
1051
                      Molcrys%Mol(k)%lS_TLS(i,j)=Molcrys%Mol(k)%lS_TLS(i,j)-1
1052
                   end if
1053
                end do
1054
             end do
1055

    
1056
             if (Molcrys%Mol(k)%natoms <=0) cycle
1057

    
1058
             do i=1,Molcrys%Mol(k)%natoms
1059
                do j=1,3
1060
                   if (MolCrys%Mol(k)%lI_coor(j,i) == N) then
1061
                      MolCrys%Mol(k)%lI_coor(j,i)=0
1062
                      MolCrys%Mol(k)%mI_coor(j,i)=0.0
1063
                      deleted=.true.
1064
                   end if
1065
                end do
1066

    
1067
                if (MolCrys%Mol(k)%lbiso(i) == N) then
1068
                   MolCrys%Mol(k)%lbiso(i)=0
1069
                   MolCrys%Mol(k)%mbiso(i)=0.0
1070
                   deleted=.true.
1071
                end if
1072

    
1073
                if (MolCrys%Mol(k)%locc(i) == N) then
1074
                   MolCrys%Mol(k)%locc(i)=0
1075
                   MolCrys%Mol(k)%mocc(i)=0.0
1076
                   deleted=.true.
1077
                end if
1078
             end do
1079

    
1080
             do i=1,Molcrys%Mol(k)%natoms
1081
                do j=1,3
1082
                   if (MolCrys%Mol(k)%lI_coor(j,i) > N) then
1083
                      MolCrys%Mol(k)%lI_coor(j,i)=MolCrys%Mol(k)%lI_coor(j,i)-1
1084
                   end if
1085
                end do
1086

    
1087
                if (MolCrys%Mol(k)%lbiso(i) > N) then
1088
                   MolCrys%Mol(k)%lbiso(i)=MolCrys%Mol(k)%lbiso(i)-1
1089
                end if
1090

    
1091
                if (MolCrys%Mol(k)%locc(i) > N) then
1092
                   MolCrys%Mol(k)%locc(i)=MolCrys%Mol(k)%locc(i)-1
1093
                end if
1094
             end do
1095

    
1096
          end do
1097
       end if
1098

    
1099
       !---- Updating V_Vectors ----!
1100
       if (deleted) then
1101
          do i=N+1,Np_refi
1102
             V_Vec(i-1)=V_Vec(i)
1103
             V_Name(i-1)=V_Name(i)
1104
             V_Bounds(:,i-1)=V_Bounds(:,i)
1105
             V_BCon(i-1)=V_BCon(i)
1106
             V_List(i-1)=V_List(i)
1107
          end do
1108
          V_Vec(np_refi)=0.0
1109
          V_Name(np_refi)=" "
1110
          V_Bounds(:,np_refi)=0.0
1111
          V_BCon(np_refi)=0
1112
          V_List(np_refi)=0
1113

    
1114
          np_refi=np_refi-1
1115
       end if
1116

    
1117
       return
1118
    End Subroutine Delete_RefCodes_MolCrys
1119

    
1120
    !!--++
1121
    !!--++ Subroutine Delete_RefCodes_Molec(N,Molec)
1122
    !!--++    integer,             intent(in)     :: N
1123
    !!--++    type(molecule_type), intent(in out) :: Molec
1124
    !!--++
1125
    !!--++ Overloaded
1126
    !!--++ Delete the number of Refinable Parameter (N) on the list
1127
    !!--++
1128
    !!--++ Update: March - 2005
1129
    !!
1130
    Subroutine Delete_RefCodes_Molec(N,Molec)
1131
       !---- Arguments ----!
1132
       integer,             intent(in)     :: N
1133
       type(molecule_type), intent(in out) :: Molec
1134

    
1135
       !---- Local Variables ----!
1136
       logical :: deleted
1137
       integer :: i,j
1138

    
1139
       deleted=.false.
1140

    
1141
       do j=1,3
1142
          if (Molec%lxcentre(j) == N) then
1143
             Molec%lxcentre(j)=0
1144
             Molec%mxcentre(j)=0.0
1145
             deleted=.true.
1146
          end if
1147
       end do
1148

    
1149
       do j=1,3
1150
          if (Molec%lOrient(j) == N) then
1151
             Molec%lOrient(j)=0
1152
             Molec%mOrient(j)=0.0
1153
             deleted=.true.
1154
          end if
1155
       end do
1156

    
1157
       do j=1,6
1158
          if (Molec%lT_TLS(j) == N) then
1159
             Molec%lT_TLS(j)=0
1160
             Molec%mT_TLS(j)=0.0
1161
             deleted=.true.
1162
          end if
1163
       end do
1164

    
1165
       do j=1,6
1166
          if (Molec%lL_TLS(j) == N) then
1167
             Molec%lL_TLS(j)=0
1168
             Molec%mL_TLS(j)=0.0
1169
             deleted=.true.
1170
          end if
1171
       end do
1172

    
1173
       do i=1,3
1174
          do j=1,3
1175
             if (Molec%lS_TLS(i,j) == N) then
1176
                Molec%lS_TLS(i,j)=0
1177
                Molec%mS_TLS(i,j)=0.0
1178
                deleted=.true.
1179
             end if
1180
          end do
1181
       end do
1182

    
1183
       !---- Updating ----!
1184
       do j=1,3
1185
          if (Molec%lxcentre(j) > N) then
1186
             Molec%lxcentre(j)=Molec%lxcentre(j)-1
1187
          end if
1188
       end do
1189

    
1190
       do j=1,3
1191
          if (Molec%lOrient(j) > N) then
1192
             Molec%lOrient(j)=Molec%lOrient(j)-1
1193
          end if
1194
       end do
1195

    
1196
       do j=1,6
1197
          if (Molec%lT_TLS(j) > N) then
1198
             Molec%lT_TLS(j)=Molec%lT_TLS(j)-1
1199
          end if
1200
       end do
1201

    
1202
       do j=1,6
1203
          if (Molec%lL_TLS(j) > N) then
1204
             Molec%lL_TLS(j)=Molec%lL_TLS(j)-1
1205
          end if
1206
       end do
1207

    
1208
       do i=1,3
1209
          do j=1,3
1210
             if (Molec%lS_TLS(i,j) > N) then
1211
                Molec%lS_TLS(i,j)=Molec%lS_TLS(i,j)-1
1212
             end if
1213
          end do
1214
       end do
1215

    
1216
       if (molec%natoms <=0) return
1217

    
1218
       do i=1,Molec%Natoms
1219
          do j=1,3
1220
             if (Molec%lI_coor(j,i) == N) then
1221
                Molec%lI_coor(j,i)=0
1222
                Molec%mI_coor(j,i)=0.0
1223
                deleted=.true.
1224
             end if
1225
          end do
1226

    
1227
          if (Molec%lbiso(i) == N) then
1228
             Molec%lbiso(i)=0
1229
             Molec%mbiso(i)=0.0
1230
             deleted=.true.
1231
          end if
1232

    
1233
          if (Molec%locc(i) == N) then
1234
             Molec%locc(i)=0
1235
             Molec%mocc(i)=0.0
1236
             deleted=.true.
1237
          end if
1238
       end do
1239

    
1240
       !---- Updating ----!
1241
       do i=1,Molec%Natoms
1242
          do j=1,3
1243
             if (Molec%lI_coor(j,i) > N) then
1244
                Molec%lI_coor(j,i)=Molec%lI_coor(j,i)-1
1245
             end if
1246
          end do
1247

    
1248
          if (Molec%lbiso(i) > N) then
1249
             Molec%lbiso(i)=Molec%lbiso(i)-1
1250
          end if
1251

    
1252
          if (Molec%locc(i) > N) then
1253
             Molec%locc(i)=Molec%locc(i)-1
1254
          end if
1255
       end do
1256

    
1257
       !---- Updating V_Vectors ----!
1258
       if (deleted) then
1259
          do i=N+1,Np_refi
1260
             V_Vec(i-1)=V_Vec(i)
1261
             V_Name(i-1)=V_Name(i)
1262
             V_Bounds(:,i-1)=V_Bounds(:,i)
1263
             V_BCon(i-1)=V_BCon(i)
1264
             V_List(i-1)=V_List(i)
1265
          end do
1266
          V_Vec(np_refi)=0.0
1267
          V_Name(np_refi)=" "
1268
          V_Bounds(:,np_refi)=0.0
1269
          V_BCon(np_refi)=0
1270
          V_List(np_refi)=0
1271

    
1272
          np_refi=np_refi-1
1273
       end if
1274

    
1275
       return
1276
    End Subroutine Delete_RefCodes_Molec
1277

    
1278
    !!--++
1279
    !!--++ Subroutine Delete_RefCodes_Magdom(N, Mag_Dom)
1280
    !!--++    integer,              intent(in)          :: N
1281
    !!--++    type(Magnetic_Domain_type),intent(in out) :: Mag_Dom
1282
    !!--++
1283
    !!--++ Delete the number of Refinable Parameters (N) on the list
1284
    !!--++ related to magnetic domains
1285
    !!--++ Created: February - 2012
1286
    !!
1287
    Subroutine Delete_RefCodes_Magdom(N, Mag_Dom)
1288
       !---- Arguments ----!
1289
       integer,              intent(in)           :: N
1290
       type(Magnetic_Domain_type), intent(in out) :: Mag_Dom
1291

    
1292
       !---- Local Variables ----!
1293
       logical :: deleted
1294
       integer :: i,j,ich
1295

    
1296
       deleted=.false.
1297

    
1298
       !---- Check is chirality is present ----!
1299
       if (Mag_Dom%chir) then
1300
        ich=2
1301
       else
1302
        ich=1        
1303
       end if
1304

    
1305
       !---- Eliminate N Parameter ----!
1306
        do i=1,Mag_Dom%nd
1307
         do j=1,ich
1308
          if (Mag_Dom%Lpop(j,i) == N) then
1309
              Mag_Dom%Mpop(j,i)=0.0
1310
              Mag_Dom%Lpop(j,i)=0
1311
              deleted=.true.
1312
          end if
1313
         end do
1314
        end do
1315

    
1316
       !---- Updating Variables ----!
1317
        do i=1,Mag_Dom%nd
1318
         do j=1,ich
1319
          if (Mag_Dom%Lpop(j,i) > N) then
1320
              Mag_Dom%Lpop(j,i)=Mag_Dom%Lpop(j,i)-1
1321
          end if
1322
         end do
1323
        end do
1324

    
1325
       !---- Updating V_Vectors ----!
1326
       if (deleted) then
1327
          do i=N+1,Np_refi
1328
             V_Vec(i-1)=V_Vec(i)
1329
             V_Name(i-1)=V_Name(i)
1330
             V_Bounds(:,i-1)=V_Bounds(:,i)
1331
             V_BCon(i-1)=V_BCon(i)
1332
             V_List(i-1)=V_List(i)
1333
          end do
1334
          V_Vec(np_refi)=0.0
1335
          V_Name(np_refi)=" "
1336
          V_Bounds(:,np_refi)=0.0
1337
          V_BCon(np_refi)=0
1338
          V_List(np_refi)=0
1339

    
1340
          np_refi=np_refi-1
1341
       end if
1342

    
1343
       return
1344
    End Subroutine Delete_RefCodes_Magdom
1345

    
1346
    !!--++
1347
    !!--++ Subroutine Fill_RefCodes(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom/MolCrys/Molec,Spg)
1348
    !!--++    integer,                       intent(in)     :: Key
1349
    !!--++    character(len=*),              intent(in)     :: Dire
1350
    !!--++    integer,                       intent(in)     :: Na
1351
    !!--++    integer,                       intent(in)     :: Nb
1352
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1353
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1354
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1355
    !!--++    integer,                       intent(in)     :: Ic
1356
    !!--++    type(Atom_List_Type),          intent(in out) :: FAtom
1357
    !!--++    or
1358
    !!--++    type(molecular_Crystal_type),  intent(in out) :: MolCrys
1359
    !!--++    or
1360
    !!--++    type(molecule_type),           intent(in out) :: Molec
1361
    !!--++    type(space_group_type),        intent(in)     :: Spg
1362
    !!--++
1363
    !!--++    (Private)
1364
    !!--++    Write on Vectors the Information for Free Atoms
1365
    !!--++
1366
    !!--++ Update: March - 2005
1367
    !!
1368

    
1369
    !!--++
1370
    !!--++ Subroutine Fill_RefCodes_FAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom,Spg)
1371
    !!--++    integer,                       intent(in)     :: Key
1372
    !!--++    character(len=*),              intent(in)     :: Dire
1373
    !!--++    integer,                       intent(in)     :: Na
1374
    !!--++    integer,                       intent(in)     :: Nb
1375
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1376
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1377
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1378
    !!--++    integer,                       intent(in)     :: Ic
1379
    !!--++    type(Atom_List_Type),          intent(in out) :: FAtom
1380
    !!--++    type(space_group_type),        intent(in)     :: Spg
1381
    !!--++
1382
    !!--++ Overloaded
1383
    !!--++ Write on Vectors the Information for Free Atoms
1384
    !!--++
1385
    !!--++ Update: March - 2005
1386
    !!
1387
    Subroutine Fill_RefCodes_FAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom,Spg)
1388
       !---- Arguments ----!
1389
       integer,                       intent(in)     :: Key
1390
       character(len=*),              intent(in)     :: Dire
1391
       integer,                       intent(in)     :: Na
1392
       integer,                       intent(in)     :: Nb
1393
       real(kind=cp),                 intent(in)     :: Xl
1394
       real(kind=cp),                 intent(in)     :: Xu
1395
       real(kind=cp),                 intent(in)     :: Xs
1396
       integer,                       intent(in)     :: Ic
1397
       type(Atom_List_Type),          intent(in out) :: FAtom
1398
       type(space_group_type),        intent(in)     :: Spg
1399

    
1400
       !---- Local variables ----!
1401
       integer           :: j,nc
1402

    
1403
       call init_err_refcodes()
1404
       if (Na <= 0) then
1405
          err_refcodes=.true.
1406
          ERR_RefCodes_Mess="Number of atom no defined"
1407
          return
1408
       end if
1409

    
1410
       select case (dire)
1411
          !---- FIX Directive ----!
1412
          case ("fix")
1413

    
1414
             select case (key)
1415
                case (0)
1416
                   !---- nb must be different zero ----!
1417
                   select case (nb)
1418
                      case (0)
1419
                         err_refcodes=.true.
1420
                         ERR_RefCodes_Mess="Option not defined"
1421
                         return
1422

    
1423
                      case ( 1:3)
1424
                         !---- X_, Y_, Z_ ----!
1425
                         if (FAtom%atom(na)%lx(nb) /=0) then
1426
                            nc=FAtom%atom(na)%lx(nb)
1427
                            call Delete_RefCodes(nc,FAtom)
1428
                         end if
1429

    
1430
                      case ( 4)
1431
                         !---- Biso_ ----!
1432
                         if (FAtom%atom(na)%lbiso /=0) then
1433
                            nc=FAtom%atom(na)%lbiso
1434
                            call Delete_RefCodes(nc,fatom)
1435
                         end if
1436

    
1437
                      case ( 5)
1438
                         !---- Occ_ ----!
1439
                         if (FAtom%atom(na)%locc /=0) then
1440
                            nc=FAtom%atom(na)%locc
1441
                            call Delete_RefCodes(nc,fatom)
1442
                         end if
1443

    
1444
                      case ( 6:11)
1445
                         !---- B11_,...,B23_ ----!
1446
                         if (FAtom%atom(na)%lu(nb-5) /=0) then
1447
                            nc=FAtom%atom(na)%lu(nb-5)
1448
                            call Delete_RefCodes(nc,fatom)
1449
                         end if
1450

    
1451
                      case (12)
1452
                         !---- Banis_ ----!
1453
                         do j=1,6
1454
                            if (FAtom%atom(na)%lu(j) /=0) then
1455
                               nc=FAtom%atom(na)%lu(j)
1456
                               call Delete_RefCodes(nc,fatom)
1457
                            end if
1458
                         end do
1459

    
1460
                      case (13:)
1461
                         err_refcodes=.true.
1462
                         ERR_RefCodes_Mess="Option not defined for this type of variable "
1463
                         return
1464
                   end select ! nb
1465

    
1466
                case (1)
1467
                   !---- XYZ ----!
1468
                   do j=1,3
1469
                      if (FAtom%atom(na)%lx(j) /=0) then
1470
                         nc=FAtom%atom(na)%lx(j)
1471
                         call Delete_RefCodes(nc,fatom)
1472
                      end if
1473
                   end do
1474

    
1475
                case (2)
1476
                   !---- OCC ----!
1477
                   if (FAtom%atom(na)%locc /=0) then
1478
                      nc=FAtom%atom(na)%locc
1479
                      call Delete_RefCodes(nc,fatom)
1480
                   end if
1481

    
1482
                case (3)
1483
                   !---- BIS ----!
1484
                   if (FAtom%atom(na)%lbiso /=0) then
1485
                      nc=FAtom%atom(na)%lbiso
1486
                      call Delete_RefCodes(nc,fatom)
1487
                   end if
1488

    
1489
                case (4)
1490
                  !---- BAN ----!
1491
                  do j=1,6
1492
                     if (FAtom%atom(na)%lu(j) /=0) then
1493
                        nc=FAtom%atom(na)%lu(j)
1494
                        call Delete_RefCodes(nc,fatom)
1495
                     end if
1496
                  end do
1497

    
1498
                case (5)
1499
                   !---- ALL ----!
1500
                   do j=1,3
1501
                      if (FAtom%atom(na)%lx(j) /=0) then
1502
                         nc=FAtom%atom(na)%lx(j)
1503
                         call Delete_RefCodes(nc,fatom)
1504
                      end if
1505
                   end do
1506
                   if (FAtom%atom(na)%locc /=0) then
1507
                      nc=FAtom%atom(na)%locc
1508
                      call Delete_RefCodes(nc,fatom)
1509
                   end if
1510
                   if (FAtom%atom(na)%lbiso /=0) then
1511
                      nc=FAtom%atom(na)%lbiso
1512
                      call Delete_RefCodes(nc,fatom)
1513
                   end if
1514
                   do j=1,6
1515
                      if (FAtom%atom(na)%lu(j) /=0) then
1516
                         nc=FAtom%atom(na)%lu(j)
1517
                         call Delete_RefCodes(nc,fatom)
1518
                      end if
1519
                   end do
1520

    
1521
                case (6:)
1522
                   err_refcodes=.true.
1523
                   ERR_RefCodes_Mess="Incompatible information for this type of variable "
1524
                   return
1525
             end select
1526

    
1527
          !---- VARY Directive ----!
1528
          case ("var")
1529

    
1530
             select case (key)
1531
                case (0)
1532

    
1533
                   !---- nb must be different zero ----!
1534
                   select case (nb)
1535
                      case (0)
1536
                         err_refcodes=.true.
1537
                         ERR_RefCodes_Mess="Option not defined"
1538
                         return
1539

    
1540
                      case ( 1:3)
1541
                         !---- X_, Y_, Z_ ----!
1542
                         if (FAtom%atom(na)%lx(nb) ==0) then
1543
                            FAtom%atom(na)%mx(nb)=1.0
1544
                            call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1545
                                                 FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1546
                            if (FAtom%atom(na)%lx(nb) == np_refi) then
1547
                               V_Vec(np_refi)=FAtom%atom(na)%x(nb)
1548
                               V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1549
                               V_Bounds(1,np_refi)=xl
1550
                               V_Bounds(2,np_refi)=xu
1551
                               V_Bounds(3,np_refi)=xs
1552
                               V_BCon(np_refi)=ic
1553
                               V_list(np_refi)=na
1554
                            else
1555
                               np_refi=np_refi-1
1556
                            end if
1557
                         end if
1558

    
1559
                      case ( 4)
1560
                         !---- Biso_ ----!
1561
                         if (FAtom%atom(na)%lbiso ==0) then
1562
                            np_refi=np_refi+1
1563
                            V_Vec(np_refi)=FAtom%atom(na)%biso
1564
                            V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1565
                            FAtom%atom(na)%mbiso=1.0
1566
                            FAtom%atom(na)%lbiso=np_refi
1567
                            V_Bounds(1,np_refi)=xl
1568
                            V_Bounds(2,np_refi)=xu
1569
                            V_Bounds(3,np_refi)=xs
1570
                            V_BCon(np_refi)=ic
1571
                            V_list(np_refi)=na
1572
                         end if
1573

    
1574
                      case ( 5)
1575
                         !---- Occ_ ----!
1576
                         if (FAtom%atom(na)%locc ==0) then
1577
                            np_refi=np_refi+1
1578
                            V_Vec(np_refi)=FAtom%atom(na)%occ
1579
                            V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1580
                            FAtom%atom(na)%mocc=1.0
1581
                            FAtom%atom(na)%locc=np_refi
1582
                            V_Bounds(1,np_refi)=xl
1583
                            V_Bounds(2,np_refi)=xu
1584
                            V_Bounds(3,np_refi)=xs
1585
                            V_BCon(np_refi)=ic
1586
                            V_list(np_refi)=na
1587
                         end if
1588

    
1589
                      case ( 6:11)
1590
                         !---- B11_,...,B23_ ----!
1591
                         if (FAtom%atom(na)%lu(nb-5) ==0) then
1592
                            FAtom%atom(na)%mu(nb-5)=1.0
1593
                            call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1594
                                                 np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1595
                            if (FAtom%atom(na)%lu(nb-5) == np_refi) then
1596
                               V_Vec(np_refi)=FAtom%atom(na)%u(nb-5)
1597
                               V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1598
                               V_Bounds(1,np_refi)=xl
1599
                               V_Bounds(2,np_refi)=xu
1600
                               V_Bounds(3,np_refi)=xs
1601
                               V_BCon(np_refi)=ic
1602
                               V_list(np_refi)=na
1603
                            else
1604
                               np_refi=np_refi-1
1605
                            end if
1606
                         end if
1607

    
1608
                      case (12)
1609
                         !---- Banis_ ----!
1610
                         do j=1,6
1611
                            if (FAtom%atom(na)%lu(j) ==0) then
1612
                               FAtom%atom(na)%mu(j)=1.0
1613
                               call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1614
                                                    np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1615
                               if (FAtom%atom(na)%lu(j) == np_refi) then
1616
                                  V_Vec(np_refi)=FAtom%atom(na)%u(j)
1617
                                  V_Name(np_refi)=trim(code_nam(5+j))//trim(FAtom%atom(na)%lab)
1618
                                  V_Bounds(1,np_refi)=xl
1619
                                  V_Bounds(2,np_refi)=xu
1620
                                  V_Bounds(3,np_refi)=xs
1621
                                  V_BCon(np_refi)=ic
1622
                                  V_list(np_refi)=na
1623
                               else
1624
                                  np_refi=np_refi-1
1625
                               end if
1626
                            end if
1627
                         end do
1628

    
1629
                      case (13:)
1630
                         err_refcodes=.true.
1631
                         ERR_RefCodes_Mess="Option Not defined by this type of variables"
1632
                         return
1633

    
1634
                   end select ! nb
1635

    
1636
                case (1)
1637
                   !---- XYZ ----!
1638
                   do j=1,3
1639
                      if (FAtom%atom(na)%lx(j) ==0) then
1640
                         FAtom%atom(na)%mx(j)=1.0
1641
                         call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1642
                                              FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1643
                         if (FAtom%atom(na)%lx(j) == np_refi) then
1644
                            V_Vec(np_refi)=FAtom%atom(na)%x(j)
1645
                            V_Name(np_refi)=trim(code_nam(j))//trim(FAtom%atom(na)%lab)
1646
                            V_Bounds(1,np_refi)=xl
1647
                            V_Bounds(2,np_refi)=xu
1648
                            V_Bounds(3,np_refi)=xs
1649
                            V_BCon(np_refi)=ic
1650
                            V_list(np_refi)=na
1651
                         else
1652
                            np_refi=np_refi-1
1653
                         end if
1654
                      end if
1655
                   end do
1656

    
1657
                case (2)
1658
                   !---- OCC ----!
1659
                   if (FAtom%atom(na)%locc ==0) then
1660
                      np_refi=np_refi+1
1661
                      V_Vec(np_refi)=FAtom%atom(na)%occ
1662
                      V_Name(np_refi)=trim(code_nam(5))//trim(FAtom%atom(na)%lab)
1663
                      FAtom%atom(na)%mocc=1.0
1664
                      FAtom%atom(na)%locc=np_refi
1665
                      V_Bounds(1,np_refi)=xl
1666
                      V_Bounds(2,np_refi)=xu
1667
                      V_Bounds(3,np_refi)=xs
1668
                      V_BCon(np_refi)=ic
1669
                      V_list(np_refi)=na
1670
                   end if
1671

    
1672
                case (3)
1673
                   !---- BIS ----!
1674
                   if (FAtom%atom(na)%lbiso ==0) then
1675
                      np_refi=np_refi+1
1676
                      V_Vec(np_refi)=FAtom%atom(na)%biso
1677
                      V_Name(np_refi)=trim(code_nam(4))//trim(FAtom%atom(na)%lab)
1678
                      FAtom%atom(na)%mbiso=1.0
1679
                      FAtom%atom(na)%lbiso=np_refi
1680
                      V_Bounds(1,np_refi)=xl
1681
                      V_Bounds(2,np_refi)=xu
1682
                      V_Bounds(3,np_refi)=xs
1683
                      V_BCon(np_refi)=ic
1684
                      V_list(np_refi)=na
1685
                   end if
1686

    
1687
                case (4)
1688
                   !---- BAN ----!
1689
                   do j=1,6
1690
                      if (FAtom%atom(na)%lu(j) ==0) then
1691
                         FAtom%atom(na)%mu(j)=1.0
1692
                         call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1693
                                              np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1694
                         if (FAtom%atom(na)%lu(j) == np_refi) then
1695
                            V_Vec(np_refi)=FAtom%atom(na)%u(j)
1696
                            V_Name(np_refi)=trim(code_nam(5+j))//trim(FAtom%atom(na)%lab)
1697
                            V_Bounds(1,np_refi)=xl
1698
                            V_Bounds(2,np_refi)=xu
1699
                            V_Bounds(3,np_refi)=xs
1700
                            V_BCon(np_refi)=ic
1701
                            V_list(np_refi)=na
1702
                         else
1703
                            np_refi=np_refi-1
1704
                         end if
1705
                      end if
1706
                   end do
1707

    
1708
                case (5)
1709
                   !---- ALL ----!
1710
                   do j=1,3
1711
                      if (FAtom%atom(na)%lx(j) ==0) then
1712
                         FAtom%atom(na)%mx(j)=1.0
1713
                         call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1714
                                              FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1715
                         if (FAtom%atom(na)%lx(j) == np_refi) then
1716
                            V_Vec(np_refi)=FAtom%atom(na)%x(j)
1717
                            V_Name(np_refi)=trim(code_nam(j))//trim(FAtom%atom(na)%lab)
1718
                            V_Bounds(1,np_refi)=xl
1719
                            V_Bounds(2,np_refi)=xu
1720
                            V_Bounds(3,np_refi)=xs
1721
                            V_BCon(np_refi)=ic
1722
                            V_list(np_refi)=na
1723
                         else
1724
                            np_refi=np_refi-1
1725
                         end if
1726
                      end if
1727
                   end do
1728
                   if (FAtom%atom(na)%locc ==0) then
1729
                      np_refi=np_refi+1
1730
                      V_Vec(np_refi)=FAtom%atom(na)%occ
1731
                      V_Name(np_refi)=trim(mcode_nam(5))//trim(FAtom%atom(na)%lab)
1732
                      FAtom%atom(na)%mocc=1.0
1733
                      FAtom%atom(na)%locc=np_refi
1734
                      V_Bounds(1,np_refi)=xl
1735
                      V_Bounds(2,np_refi)=xu
1736
                      V_Bounds(3,np_refi)=xs
1737
                      V_BCon(np_refi)=ic
1738
                      V_list(np_refi)=na
1739
                   end if
1740
                   if (FAtom%atom(na)%lbiso ==0) then
1741
                      np_refi=np_refi+1
1742
                      V_Vec(np_refi)=FAtom%atom(na)%biso
1743
                      V_Name(np_refi)=trim(mcode_nam(4))//trim(FAtom%atom(na)%lab)
1744
                      FAtom%atom(na)%mbiso=1.0
1745
                      FAtom%atom(na)%lbiso=np_refi
1746
                      V_Bounds(1,np_refi)=xl
1747
                      V_Bounds(2,np_refi)=xu
1748
                      V_Bounds(3,np_refi)=xs
1749
                      V_BCon(np_refi)=ic
1750
                      V_list(np_refi)=na
1751
                   end if
1752
                   do j=1,6
1753
                      if (FAtom%atom(na)%lu(j) ==0) then
1754
                         FAtom%atom(na)%mu(j)=1.0
1755
                         call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1756
                                              np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1757
                         if (FAtom%atom(na)%lu(j) == np_refi) then
1758
                            V_Vec(np_refi)=FAtom%atom(na)%u(j)
1759
                            V_Name(np_refi)=trim(mcode_nam(5+j))//trim(FAtom%atom(na)%lab)
1760
                            V_Bounds(1,np_refi)=xl
1761
                            V_Bounds(2,np_refi)=xu
1762
                            V_Bounds(3,np_refi)=xs
1763
                            V_BCon(np_refi)=ic
1764
                            V_list(np_refi)=na
1765
                         else
1766
                            np_refi=np_refi-1
1767
                         end if
1768
                      end if
1769
                   end do
1770

    
1771
                case(6:)
1772
                   err_refcodes=.true.
1773
                   ERR_RefCodes_Mess="Option Not defined by this type of variables"
1774
                   return
1775
             end select
1776
       end select
1777

    
1778
       return
1779
    End Subroutine Fill_RefCodes_FAtom
1780

    
1781
    !!--++
1782
    !!--++ Subroutine Fill_RefCodes_FmAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FmAtom,Spg)
1783
    !!--++    integer,                       intent(in)     :: Key
1784
    !!--++    character(len=*),              intent(in)     :: Dire
1785
    !!--++    integer,                       intent(in)     :: Na
1786
    !!--++    integer,                       intent(in)     :: Nb
1787
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1788
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1789
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1790
    !!--++    integer,                       intent(in)     :: Ic
1791
    !!--++    type(mAtom_List_Type),         intent(in out) :: FmAtom
1792
    !!--++    type(space_group_type),        intent(in)     :: Spg
1793
    !!--++
1794
    !!--++ Write on Vectors the Information for Free Magnetic Atoms
1795
    !!--++ magnetic clone of Fill_RefCodes_FAtom
1796
    !!--++ Created: December - 2011
1797
    !!--++
1798
    Subroutine Fill_RefCodes_FmAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FmAtom,Spg)
1799
       !---- Arguments ----!
1800
       integer,                       intent(in)     :: Key
1801
       character(len=*),              intent(in)     :: Dire
1802
       integer,                       intent(in)     :: Na
1803
       integer,                       intent(in)     :: Nb
1804
       real(kind=cp),                 intent(in)     :: Xl
1805
       real(kind=cp),                 intent(in)     :: Xu
1806
       real(kind=cp),                 intent(in)     :: Xs
1807
       integer,                       intent(in)     :: Ic
1808
       type(mAtom_List_Type),         intent(in out) :: FmAtom
1809
       type(space_group_type),        intent(in)     :: Spg
1810

    
1811
       !---- Local variables ----!
1812
       integer           :: j,nc,ik
1813

    
1814
       call init_err_refcodes()
1815
       if (Na <= 0) then
1816
          err_refcodes=.true.
1817
          ERR_RefCodes_Mess="Number of atom no defined"
1818
          return
1819
       end if
1820

    
1821
       !---- Get im, number of the magnetic matrices/irrep set
1822
       ik=FmAtom%atom(na)%nvk
1823

    
1824
       select case (dire)
1825
          !---- FIX Directive ----!
1826
          case ("fix")
1827

    
1828
             select case (key)
1829
                case (0)
1830
                   !---- nb must be different zero ----!
1831
                   select case (nb)
1832
                      case (0)
1833
                         err_refcodes=.true.
1834
                         ERR_RefCodes_Mess="Option not defined"
1835
                         return
1836

    
1837
                      case ( 1:3)
1838
                         !---- Rx_, Ry_, Rz_ ----!
1839
                         if (FmAtom%atom(na)%lSkR(nb,ik) /=0) then
1840
                            nc=FmAtom%atom(na)%lSkR(nb,ik)
1841
                            call Delete_RefCodes(nc,FmAtom)
1842
                         end if
1843

    
1844
                      case ( 4:6)
1845
                         !---- Ix_, Iy_, Iz_ ----!
1846
                         if (FmAtom%atom(na)%lSkI(nb-3,ik) /=0) then
1847
                            nc=FmAtom%atom(na)%lSkI(nb-3,ik)
1848
                            call Delete_RefCodes(nc,FmAtom)
1849
                         end if
1850

    
1851
                      case ( 7:9)
1852
                         !---- Rm_, Rphi_, Rth_ ----!
1853
                         if (FmAtom%atom(na)%lSkR(nb-6,ik) /=0) then
1854
                            nc=FmAtom%atom(na)%lSkR(nb-6,ik)
1855
                            call Delete_RefCodes(nc,FmAtom)
1856
                         end if
1857

    
1858
                      case (10:12)
1859
                         !---- Im_, Iphi_, Ith_ ----!
1860
                         if (FmAtom%atom(na)%lSkI(nb-9,ik) /=0) then
1861
                            nc=FmAtom%atom(na)%lSkI(nb-9,ik)
1862
                            call Delete_RefCodes(nc,FmAtom)
1863
                         end if
1864

    
1865
                      case (13)
1866
                         !---- MagPh_ ----!
1867
                            if (FmAtom%atom(na)%lmphas(ik) /=0) then
1868
                               nc=FmAtom%atom(na)%lmphas(ik)
1869
                               call Delete_RefCodes(nc,FmAtom)
1870
                            end if
1871

    
1872
                      case (14:22)
1873
                         !---- C1_,..C12_ ----!
1874
                            if (FmAtom%atom(na)%lbas(nb-13,ik) /=0) then
1875
                               nc=FmAtom%atom(na)%lbas(nb-13,ik)
1876
                               call Delete_RefCodes(nc,FmAtom)
1877
                            end if
1878

    
1879
                      case (23)
1880
                         err_refcodes=.true.
1881
                         ERR_RefCodes_Mess="Option not defined for this type of variable "
1882
                         return
1883
                   end select ! nb
1884

    
1885
                case (1)
1886
                   !---- Rxyz ----!
1887
                   do j=1,3
1888
                      if (FmAtom%atom(na)%lSkR(j,ik) /=0) then
1889
                         nc=FmAtom%atom(na)%lSkR(j,ik)
1890
                         call Delete_RefCodes(nc,FmAtom)
1891
                      end if
1892
                   end do
1893

    
1894
                case (2)
1895
                   !---- Ixyz ----!
1896
                   do j=1,3
1897
                      if (FmAtom%atom(na)%lSkI(j,ik) /=0) then
1898
                         nc=FmAtom%atom(na)%lSkI(j,ik)
1899
                         call Delete_RefCodes(nc,FmAtom)
1900
                      end if
1901
                   end do
1902

    
1903
                case (3)
1904
                   !---- Mxyz ----!
1905
                   do j=1,3
1906
                      if (FmAtom%atom(na)%lSkR(j,ik) /=0) then
1907
                         nc=FmAtom%atom(na)%lSkR(j,ik)
1908
                         call Delete_RefCodes(nc,FmAtom)
1909
                      end if
1910
                   end do
1911
                   do j=1,3
1912
                      if (FmAtom%atom(na)%lSkI(j,ik) /=0) then
1913
                         nc=FmAtom%atom(na)%lSkI(j,ik)
1914
                         call Delete_RefCodes(nc,FmAtom)
1915
                      end if
1916
                   end do
1917

    
1918
                case (4:)
1919
                   err_refcodes=.true.
1920
                   ERR_RefCodes_Mess="Incompatible information for this type of variable "
1921
                   return
1922
             end select
1923

    
1924
          !---- VARY Directive ----!
1925
          case ("var")
1926

    
1927
              select case (key)
1928
                case (0)
1929

    
1930
                   !---- nb must be different zero ----!
1931
                   select case (nb)
1932
                      case (0)
1933
                         err_refcodes=.true.
1934
                         ERR_RefCodes_Mess="Option not defined"
1935
                         return
1936

    
1937
                      case ( 1:3)
1938
                         !---- Rx_, Ry_, Rz_ ----!
1939
                         if (FmAtom%atom(na)%lSkR(nb,ik) ==0) then
1940

    
1941
                            np_refi=np_refi+1
1942
                            FmAtom%atom(na)%lSkR(nb,ik)=np_refi
1943
                            FmAtom%atom(na)%mSkR(nb,ik)=1.0
1944

    
1945
                            if (FmAtom%atom(na)%lSkR(nb,ik) == np_refi) then
1946
                               V_Vec(np_refi)=FmAtom%atom(na)%SkR(nb,ik)
1947
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1948
                               V_Bounds(1,np_refi)=xl
1949
                               V_Bounds(2,np_refi)=xu
1950
                               V_Bounds(3,np_refi)=xs
1951
                               V_BCon(np_refi)=ic
1952
                               V_list(np_refi)=na
1953
                            else
1954
                               np_refi=np_refi-1
1955
                            end if
1956
                         end if
1957

    
1958
                      case ( 4:6)
1959
                         !---- Ix_, Iy_, Iz_ ----!
1960
                         if (FmAtom%atom(na)%lSkI(nb-3,ik) ==0) then
1961

    
1962
                            np_refi=np_refi+1
1963
                            FmAtom%atom(na)%lSkI(nb-3,ik)=np_refi
1964
                            FmAtom%atom(na)%mSkI(nb-3,ik)=1.0
1965

    
1966
                            if (FmAtom%atom(na)%lSkI(nb-3,ik) == np_refi) then
1967
                               V_Vec(np_refi)=FmAtom%atom(na)%SkI(nb-3,ik)
1968
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1969
                               V_Bounds(1,np_refi)=xl
1970
                               V_Bounds(2,np_refi)=xu
1971
                               V_Bounds(3,np_refi)=xs
1972
                               V_BCon(np_refi)=ic
1973
                               V_list(np_refi)=na
1974
                             else
1975
                               np_refi=np_refi-1
1976
                            end if
1977
                         end if
1978

    
1979
                      case ( 7:9)
1980
                         !---- Rm_, Rphi_, Rth_ ----!
1981
                         if (FmAtom%atom(na)%lSkR(nb-6,ik) ==0) then
1982

    
1983
                            np_refi=np_refi+1
1984
                            FmAtom%atom(na)%lSkR(nb-6,ik)=np_refi
1985
                            FmAtom%atom(na)%mSkR(nb-6,ik)=1.0
1986

    
1987
                            if (FmAtom%atom(na)%lSkR(nb-6,ik) == np_refi) then
1988
                               V_Vec(np_refi)=FmAtom%atom(na)%Spher_SkR(nb-6,ik)
1989
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1990
                               V_Bounds(1,np_refi)=xl
1991
                               V_Bounds(2,np_refi)=xu
1992
                               V_Bounds(3,np_refi)=xs
1993
                               V_BCon(np_refi)=ic
1994
                               V_list(np_refi)=na
1995
                            else
1996
                               np_refi=np_refi-1
1997
                            end if
1998
                         end if
1999

    
2000
                      case (10:12)
2001
                         !---- Im_, Iphi_, Ith_ ----!
2002
                         if (FmAtom%atom(na)%lSkI(nb-9,ik) ==0) then
2003

    
2004
                            np_refi=np_refi+1
2005
                            FmAtom%atom(na)%lSkI(nb-9,ik)=np_refi
2006
                            FmAtom%atom(na)%mSkI(nb-9,ik)=1.0
2007

    
2008
                            if (FmAtom%atom(na)%lSkI(nb-9,ik) == np_refi) then
2009
                               V_Vec(np_refi)=FmAtom%atom(na)%Spher_SkI(nb-9,ik)
2010
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
2011
                               V_Bounds(1,np_refi)=xl
2012
                               V_Bounds(2,np_refi)=xu
2013
                               V_Bounds(3,np_refi)=xs
2014
                               V_BCon(np_refi)=ic
2015
                               V_list(np_refi)=na
2016
                             else
2017
                               np_refi=np_refi-1
2018
                            end if
2019
                         end if
2020

    
2021
                      case (13)
2022
                         !---- MagPh_ ----!
2023
                         if (FmAtom%atom(na)%lmphas(ik) ==0) then
2024

    
2025
                            np_refi=np_refi+1
2026
                            FmAtom%atom(na)%mmphas(ik)=1.0
2027
                            FmAtom%atom(na)%lmphas(ik)=np_refi
2028

    
2029
                            V_Vec(np_refi)=FmAtom%atom(na)%mphas(ik)
2030
                            V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
2031
                            V_Bounds(1,np_refi)=xl
2032
                            V_Bounds(2,np_refi)=xu
2033
                            V_Bounds(3,np_refi)=xs
2034
                            V_BCon(np_refi)=ic
2035
                            V_list(np_refi)=na
2036
                         end if
2037

    
2038
                      case (14:22)
2039
                         !---- C1_,..C12_ ----!
2040
                         if (FmAtom%atom(na)%lbas(nb-13,ik) ==0) then
2041

    
2042
                            np_refi=np_refi+1
2043
                            FmAtom%atom(na)%lbas(nb-13,ik)=np_refi
2044
                            FmAtom%atom(na)%mbas(nb-13,ik)=1.0
2045

    
2046
                            if (FmAtom%atom(na)%lbas(nb-13,ik) == np_refi) then
2047
                               V_Vec(np_refi)=FmAtom%atom(na)%cbas(nb-13,ik)
2048
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
2049
                               V_Bounds(1,np_refi)=xl
2050
                               V_Bounds(2,np_refi)=xu
2051
                               V_Bounds(3,np_refi)=xs
2052
                               V_BCon(np_refi)=ic
2053
                               V_list(np_refi)=na
2054
                             else
2055
                               np_refi=np_refi-1
2056
                            end if
2057
                         end if
2058

    
2059
                      case (23:)
2060
                         err_refcodes=.true.
2061
                         ERR_RefCodes_Mess="Option Not defined by this type of variables"
2062
                         return
2063

    
2064
                   end select ! nb
2065

    
2066
                case (1)
2067
                   !---- Rxyz ----!
2068
                   do j=1,3
2069
                      if (FmAtom%atom(na)%lSkR(j,ik) ==0) then
2070

    
2071
                            np_refi=np_refi+1
2072
                            FmAtom%atom(na)%lSkR(j,ik)=np_refi
2073
                            FmAtom%atom(na)%mSkR(j,ik)=1.0
2074

    
2075
                         if (FmAtom%atom(na)%lSkR(j,ik) == np_refi) then
2076
                            V_Vec(np_refi)=FmAtom%atom(na)%SkR(j,ik)
2077
                            V_Name(np_refi)=trim(mcode_nam(j))//trim(FmAtom%atom(na)%lab)
2078
                            V_Bounds(1,np_refi)=xl
2079
                            V_Bounds(2,np_refi)=xu
2080
                            V_Bounds(3,np_refi)=xs
2081
                            V_BCon(np_refi)=ic
2082
                            V_list(np_refi)=na
2083
                         else
2084
                            np_refi=np_refi-1
2085
                         end if
2086
                      end if
2087
                   end do
2088

    
2089
                case (2)
2090
                   !---- Ixyz ----!
2091
                   do j=1,3
2092
                      if (FmAtom%atom(na)%lSkI(j,ik) ==0) then
2093

    
2094
                            np_refi=np_refi+1
2095
                            FmAtom%atom(na)%lSkI(j,ik)=np_refi
2096
                            FmAtom%atom(na)%mSkI(j,ik)=1.0
2097

    
2098
                         if (FmAtom%atom(na)%lSkI(j,ik) == np_refi) then
2099
                            V_Vec(np_refi)=FmAtom%atom(na)%SkI(j,ik)
2100
                            V_Name(np_refi)=trim(mcode_nam(j+3))//trim(FmAtom%atom(na)%lab)
2101
                            V_Bounds(1,np_refi)=xl
2102
                            V_Bounds(2,np_refi)=xu
2103
                            V_Bounds(3,np_refi)=xs
2104
                            V_BCon(np_refi)=ic
2105
                            V_list(np_refi)=na
2106
                         else
2107
                            np_refi=np_refi-1
2108
                         end if
2109
                      end if
2110
                   end do
2111

    
2112
                case (3)
2113
                    !---- Mxyz ----!
2114
                   do j=1,3
2115
                      if (FmAtom%atom(na)%lSkR(j,ik) ==0) then
2116

    
2117
                            np_refi=np_refi+1
2118
                            FmAtom%atom(na)%lSkR(j,ik)=np_refi
2119
                            FmAtom%atom(na)%mSkR(j,ik)=1.0
2120

    
2121
                         if (FmAtom%atom(na)%lSkR(j,ik) == np_refi) then
2122
                            V_Vec(np_refi)=FmAtom%atom(na)%SkR(j,ik)
2123
                            V_Name(np_refi)=trim(mcode_nam(j))//trim(FmAtom%atom(na)%lab)
2124
                            V_Bounds(1,np_refi)=xl
2125
                            V_Bounds(2,np_refi)=xu
2126
                            V_Bounds(3,np_refi)=xs
2127
                            V_BCon(np_refi)=ic
2128
                            V_list(np_refi)=na
2129
                         else
2130
                            np_refi=np_refi-1
2131
                         end if
2132
                      end if
2133
                   end do
2134
                   
2135
                   do j=1,3
2136
                      if (FmAtom%atom(na)%lSkI(j,ik) ==0) then
2137

    
2138
                            np_refi=np_refi+1
2139
                            FmAtom%atom(na)%lSkI(j,ik)=np_refi
2140
                            FmAtom%atom(na)%mSkI(j,ik)=1.0
2141

    
2142
                         if (FmAtom%atom(na)%lSkI(j,ik) == np_refi) then
2143
                            V_Vec(np_refi)=FmAtom%atom(na)%SkI(j,ik)
2144
                            V_Name(np_refi)=trim(mcode_nam(j+3))//trim(FmAtom%atom(na)%lab)
2145
                            V_Bounds(1,np_refi)=xl
2146
                            V_Bounds(2,np_refi)=xu
2147
                            V_Bounds(3,np_refi)=xs
2148
                            V_BCon(np_refi)=ic
2149
                            V_list(np_refi)=na
2150
                         else
2151
                            np_refi=np_refi-1
2152
                         end if
2153
                      end if
2154
                   end do
2155
 
2156
                case (4:)
2157
                   err_refcodes=.true.
2158
                   ERR_RefCodes_Mess="Option Not defined by this type of variables"
2159
                   return
2160
             end select
2161
       end select
2162

    
2163
       return
2164
    End Subroutine Fill_RefCodes_FmAtom
2165

    
2166
    !!--++
2167
    !!--++ Subroutine Fill_RefCodes_Molcrys(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molcrys,NMol)
2168
    !!--++    integer,                      intent(in)     :: Key
2169
    !!--++    character(len=*),             intent(in)     :: Dire
2170
    !!--++    integer,                      intent(in)     :: Na
2171
    !!--++    integer,                      intent(in)     :: Nb
2172
    !!--++    real(kind=cp),                intent(in)     :: Xl
2173
    !!--++    real(kind=cp),                intent(in)     :: Xu
2174
    !!--++    real(kind=cp),                intent(in)     :: Xs
2175
    !!--++    integer,                      intent(in)     :: Ic
2176
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
2177
    !!--++    integer,                      intent(in)     :: NMol
2178
    !!--++
2179
    !!--++ Overloaded
2180
    !!--++ Write on Vectors the Information for Free Atoms
2181
    !!--++
2182
    !!--++ Update: March - 2005
2183
    !!
2184
    Subroutine Fill_RefCodes_Molcrys(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molcrys,Nmol)
2185
       !---- Arguments ----!
2186
       integer,                      intent(in)     :: Key
2187
       character(len=*),             intent(in)     :: Dire
2188
       integer,                      intent(in)     :: Na
2189
       integer,                      intent(in)     :: Nb
2190
       real(kind=cp),                intent(in)     :: Xl
2191
       real(kind=cp),                intent(in)     :: Xu
2192
       real(kind=cp),                intent(in)     :: Xs
2193
       integer,                      intent(in)     :: Ic
2194
       type(molecular_Crystal_type), intent(in out) :: MolCrys
2195
       integer,                      intent(in)     :: NMol
2196

    
2197
       !---- Local variables ----!
2198
       character(len=2) :: car
2199
       integer          :: i, j, nc, naa
2200

    
2201
       call init_err_refcodes()
2202
       if (Na <= 0) then
2203
          err_refcodes=.true.
2204
          ERR_RefCodes_Mess="Number of atom no defined"
2205
          return
2206
       end if
2207

    
2208
       select case (dire)
2209
          !---- FIX Directive ----!
2210
          case ("fix")
2211

    
2212
             select case (key)
2213
                case (0)
2214
                   !---- Nb must be different zero ----!
2215
                   select case (nb)
2216
                      case (0)
2217
                         err_refcodes=.true.
2218
                         ERR_RefCodes_Mess="Option Not defined"
2219
                         return
2220

    
2221
                      case ( 1:3)
2222
                         !---- X_, Y_, Z_ ----!
2223
                         if (na <= molcrys%n_free) then
2224
                            if (molcrys%atm(na)%lx(nb) /=0) then
2225
                               nc=molcrys%atm(na)%lx(nb)
2226
                               call Delete_RefCodes(nc,MolCrys)
2227
                            end if
2228
                         else
2229
                            naa=na-molcrys%n_free
2230
                            do i=1,molcrys%n_mol
2231
                               if (naa > molcrys%mol(i)%natoms) then
2232
                                  naa=naa-molcrys%mol(i)%natoms
2233
                                  cycle
2234
                               end if
2235
                               if (molcrys%mol(i)%lI_coor(nb,naa) /=0) then
2236
                                  nc=molcrys%mol(i)%lI_coor(nb,naa)
2237
                                  call Delete_RefCodes(nc,MolCrys)
2238
                               end if
2239
                            end do
2240
                         end if
2241

    
2242
                      case ( 4)
2243
                         !---- Biso_ ----!
2244
                         if (na <= molcrys%n_free) then
2245
                            if (molcrys%atm(na)%lbiso /=0) then
2246
                               nc=molcrys%atm(na)%lbiso
2247
                               call Delete_RefCodes(nc,MolCrys)
2248
                            end if
2249
                         else
2250
                            naa=na-molcrys%n_free
2251
                            do i=1,molcrys%n_mol
2252
                               if (naa > molcrys%mol(i)%natoms) then
2253
                                  naa=naa-molcrys%mol(i)%natoms
2254
                                  cycle
2255
                               end if
2256
                               if (molcrys%mol(i)%lbiso(naa) /=0) then
2257
                                  nc=molcrys%mol(i)%lbiso(naa)
2258
                                  call Delete_RefCodes(nc,MolCrys)
2259
                               end if
2260
                            end do
2261
                         end if
2262

    
2263
                      case ( 5)
2264
                         !---- Occ_ ----!
2265
                         if (na <= molcrys%n_free) then
2266
                            if (molcrys%atm(na)%locc /=0) then
2267
                               nc=molcrys%atm(na)%locc
2268
                               call Delete_RefCodes(nc,MolCrys)
2269
                            end if
2270
                         else
2271
                            naa=na-molcrys%n_free
2272
                            do i=1,molcrys%n_mol
2273
                               if (naa > molcrys%mol(i)%natoms) then
2274
                                  naa=naa-molcrys%mol(i)%natoms
2275
                                  cycle
2276
                               end if
2277
                               if (molcrys%mol(i)%locc(naa) /=0) then
2278
                                  nc=molcrys%mol(i)%locc(naa)
2279
                                  call Delete_RefCodes(nc,MolCrys)
2280
                               end if
2281
                            end do
2282
                         end if
2283

    
2284
                      case ( 6:11)
2285
                         !---- B11_, ..., B23_ ----!
2286
                         if (na <= molcrys%n_free) then
2287
                            if (molcrys%atm(na)%lu(nb-5) /=0) then
2288
                               nc=molcrys%atm(na)%lu(nb-5)
2289
                               call Delete_RefCodes(nc,MolCrys)
2290
                            end if
2291
                         else
2292
                            err_refcodes=.true.
2293
                            ERR_RefCodes_Mess="Option Not defined"
2294
                            return
2295
                         end if
2296

    
2297
                      case (12)
2298
                         !---- Banis_ ----!
2299
                         if (na <= molcrys%n_free) then
2300
                            do j=1,6
2301
                               if (molcrys%atm(na)%lu(j) /=0) then
2302
                                  nc=molcrys%atm(na)%lu(j)
2303
                                  call Delete_RefCodes(nc,MolCrys)
2304
                               end if
2305
                            end do
2306
                         else
2307
                            err_refcodes=.true.
2308
                            ERR_RefCodes_Mess="Option Not defined"
2309
                            return
2310
                         end if
2311

    
2312
                      case (13:15)
2313
                         !---- Xc_, Yc_, Zc_ ----!
2314
                         select case (nmol)
2315
                            case (-1)
2316
                               err_refcodes=.true.
2317
                               ERR_RefCodes_Mess="Option Not defined"
2318
                               return
2319

    
2320
                            case (0)
2321
                               do i=1,molcrys%n_mol
2322
                                  if (molcrys%mol(i)%lxcentre(nb-12) /=0) then
2323
                                     nc=molcrys%mol(i)%lxcentre(nb-12)
2324
                                     call Delete_RefCodes(nc,MolCrys)
2325
                                  end if
2326
                               end do
2327

    
2328
                            case (1:)
2329
                               if (molcrys%mol(nmol)%lxcentre(nb-12) /=0) then
2330
                                  nc=molcrys%mol(nmol)%lxcentre(nb-12)
2331
                                  call Delete_RefCodes(nc,MolCrys)
2332
                               end if
2333
                         end select
2334

    
2335
                      case (16:18)
2336
                         !---- Theta_, Phi_, Chi_ ----!
2337
                         select case (nmol)
2338
                            case (-1)
2339
                               err_refcodes=.true.
2340
                               ERR_RefCodes_Mess="Option Not defined"
2341
                               return
2342

    
2343
                            case (0)
2344
                               do i=1,molcrys%n_mol
2345
                                  if (molcrys%mol(i)%lOrient(nb-15) /=0) then
2346
                                     nc=molcrys%mol(i)%lOrient(nb-15)
2347
                                     call Delete_RefCodes(nc,MolCrys)
2348
                                  end if
2349
                               end do
2350

    
2351
                            case (1:)
2352
                               if (molcrys%mol(nmol)%lOrient(nb-15) /=0) then
2353
                                  nc=molcrys%mol(nmol)%lOrient(nb-15)
2354
                                  call Delete_RefCodes(nc,MolCrys)
2355
                               end if
2356
                         end select
2357

    
2358
                      case (19:21)
2359
                         !!! Not yet implemented !!!
2360

    
2361
                   end select ! nb
2362

    
2363
                case (1)
2364
                   !---- XYZ ----!
2365
                   if (na <= molcrys%n_free) then
2366
                      do j=1,3
2367
                         if (molcrys%atm(na)%lx(j) /=0) then
2368
                            nc=molcrys%atm(na)%lx(j)
2369
                            call Delete_RefCodes(nc,MolCrys)
2370
                         end if
2371
                      end do
2372
                   else
2373
                      naa=na-molcrys%n_free
2374
                      do i=1,molcrys%n_mol
2375
                         if (naa > molcrys%mol(i)%natoms) then
2376
                            naa=naa-molcrys%mol(i)%natoms
2377
                            cycle
2378
                         end if
2379
                         do j=1,3
2380
                            if (molcrys%mol(i)%lI_coor(j,naa) /=0) then
2381
                               nc=molcrys%mol(i)%lI_coor(j,naa)
2382
                               call Delete_RefCodes(nc,MolCrys)
2383
                            end if
2384
                         end do
2385
                      end do
2386
                   end if
2387

    
2388
                case (2)
2389
                   !---- OCC ----!
2390
                   if (na <= molcrys%n_free) then
2391
                      if (molcrys%atm(na)%locc /=0) then
2392
                         nc=molcrys%atm(na)%locc
2393
                         call Delete_RefCodes(nc,MolCrys)
2394
                      end if
2395
                   else
2396
                      naa=na-molcrys%n_free
2397
                      do i=1,molcrys%n_mol
2398
                         if (naa > molcrys%mol(i)%natoms) then
2399
                            naa=naa-molcrys%mol(i)%natoms
2400
                            cycle
2401
                         end if
2402
                         if (molcrys%mol(i)%locc(naa) /=0) then
2403
                            nc=molcrys%mol(i)%locc(naa)
2404
                            call Delete_RefCodes(nc,MolCrys)
2405
                         end if
2406
                      end do
2407
                   end if
2408

    
2409
                case (3)
2410
                   !---- BIS ----!
2411
                   if (na <= molcrys%n_free) then
2412
                      if (molcrys%atm(na)%lbiso /=0) then
2413
                         nc=molcrys%atm(na)%lbiso
2414
                         call Delete_RefCodes(nc,MolCrys)
2415
                      end if
2416
                   else
2417
                      naa=na-molcrys%n_free
2418
                      do i=1,molcrys%n_mol
2419
                         if (naa > molcrys%mol(i)%natoms) then
2420
                            naa=naa-molcrys%mol(i)%natoms
2421
                            cycle
2422
                         end if
2423
                         if (molcrys%mol(i)%lbiso(naa) /=0) then
2424
                            nc=molcrys%mol(i)%lbiso(naa)
2425
                            call Delete_RefCodes(nc,MolCrys)
2426
                         end if
2427
                      end do
2428
                   end if
2429

    
2430
                case (4)
2431
                   !---- BAN ----!
2432
                   if (na <= molcrys%n_free) then
2433
                      do j=1,6
2434
                         if (molcrys%atm(na)%lu(j) /=0) then
2435
                            nc=molcrys%atm(na)%lu(j)
2436
                            call Delete_RefCodes(nc,MolCrys)
2437
                         end if
2438
                      end do
2439
                   else
2440
                      err_refcodes=.true.
2441
                      ERR_RefCodes_Mess="Option Not defined"
2442
                      return
2443
                   end if
2444

    
2445
                case (5)
2446
                   !---- ALL ----!
2447
                   if (na <= molcrys%n_free) then
2448
                      do j=1,3
2449
                         if (molcrys%atm(na)%lx(j) /=0) then
2450
                            nc=molcrys%atm(na)%lx(j)
2451
                            call Delete_RefCodes(nc,MolCrys)
2452
                         end if
2453
                      end do
2454
                   else
2455
                      naa=na-molcrys%n_free
2456
                      do i=1,molcrys%n_mol
2457
                         if (naa > molcrys%mol(i)%natoms) then
2458
                            naa=naa-molcrys%mol(i)%natoms
2459
                            cycle
2460
                         end if
2461
                         do j=1,3
2462
                            if (molcrys%mol(i)%lI_coor(j,naa) /=0) then
2463
                               nc=molcrys%mol(i)%lI_coor(j,naa)
2464
                               call Delete_RefCodes(nc,MolCrys)
2465
                            end if
2466
                         end do
2467
                      end do
2468
                   end if
2469

    
2470
                   if (na <= molcrys%n_free) then
2471
                      if (molcrys%atm(na)%locc /=0) then
2472
                         nc=molcrys%atm(na)%locc
2473
                         call Delete_RefCodes(nc,MolCrys)
2474
                      end if
2475
                   else
2476
                      naa=na-molcrys%n_free
2477
                      do i=1,molcrys%n_mol
2478
                         if (naa > molcrys%mol(i)%natoms) then
2479
                            naa=naa-molcrys%mol(i)%natoms
2480
                            cycle
2481
                         end if
2482
                         if (molcrys%mol(i)%locc(naa) /=0) then
2483
                            nc=molcrys%mol(i)%locc(naa)
2484
                            call Delete_RefCodes(nc,MolCrys)
2485
                         end if
2486
                      end do
2487
                   end if
2488

    
2489
                   if (na <= molcrys%n_free) then
2490
                      if (molcrys%atm(na)%lbiso /=0) then
2491
                         nc=molcrys%atm(na)%lbiso
2492
                         call Delete_RefCodes(nc,MolCrys)
2493
                      end if
2494
                   else
2495
                      naa=na-molcrys%n_free
2496
                      do i=1,molcrys%n_mol
2497
                         if (naa > molcrys%mol(i)%natoms) then
2498
                            naa=naa-molcrys%mol(i)%natoms
2499
                            cycle
2500
                         end if
2501
                         if (molcrys%mol(i)%lbiso(naa) /=0) then
2502
                            nc=molcrys%mol(i)%lbiso(naa)
2503
                            call Delete_RefCodes(nc,MolCrys)
2504
                         end if
2505
                      end do
2506
                   end if
2507

    
2508
                   if (na <= molcrys%n_free) then
2509
                      do j=1,6
2510
                         if (molcrys%atm(na)%lu(j) /=0) then
2511
                            nc=molcrys%atm(na)%lu(j)
2512
                            call Delete_RefCodes(nc,MolCrys)
2513
                         end if
2514
                      end do
2515
                   end if
2516

    
2517
                   select case (nmol)
2518
                      case (-1)
2519
                         err_refcodes=.true.
2520
                         ERR_RefCodes_Mess="Option Not defined"
2521
                         return
2522

    
2523
                      case (0)
2524
                         do i=1,molcrys%n_mol
2525
                            do j=1,3
2526
                               if (molcrys%mol(i)%lxcentre(j) /=0) then
2527
                                  nc=molcrys%mol(i)%lxcentre(j)
2528
                                  call Delete_RefCodes(nc,molcrys)
2529
                               end if
2530
                            end do
2531

    
2532
                            do j=1,3
2533
                               if (molcrys%mol(i)%lOrient(j) /=0) then
2534
                                  nc=molcrys%mol(i)%lOrient(j)
2535
                                  call Delete_RefCodes(nc,molcrys)
2536
                               end if
2537
                            end do
2538
                         end do
2539

    
2540
                      case (1:)
2541
                         do j=1,3
2542
                            if (molcrys%mol(nmol)%lxcentre(j) /=0) then
2543
                               nc=molcrys%mol(nmol)%lxcentre(j)
2544
                               call Delete_RefCodes(nc,molcrys)
2545
                            end if
2546
                         end do
2547

    
2548
                         do j=1,3
2549
                            if (molcrys%mol(nmol)%lOrient(j) /=0) then
2550
                               nc=molcrys%mol(nmol)%lOrient(j)
2551
                               call Delete_RefCodes(nc,molcrys)
2552
                            end if
2553
                         end do
2554
                   end select
2555

    
2556
                   !!! Not yet Implemented !!!
2557

    
2558
                case (6)
2559
                   !---- CEN ----!
2560
                   select case (nmol)
2561
                      case (-1)
2562
                         err_refcodes=.true.
2563
                         ERR_RefCodes_Mess="Option Not defined"
2564
                         return
2565

    
2566
                      case (0)
2567
                         do i=1,molcrys%n_mol
2568
                            do j=1,3
2569
                               if (molcrys%mol(i)%lxcentre(j) /=0) then
2570
                                  nc=molcrys%mol(i)%lxcentre(j)
2571
                                  call Delete_RefCodes(nc,molcrys)
2572
                               end if
2573
                            end do
2574
                         end do
2575

    
2576
                      case (1:)
2577
                         do j=1,3
2578
                            if (molcrys%mol(nmol)%lxcentre(j) /=0) then
2579
                               nc=molcrys%mol(nmol)%lxcentre(j)
2580
                               call Delete_RefCodes(nc,molcrys)
2581
                            end if
2582
                         end do
2583
                   end select
2584

    
2585
                case (7)
2586
                   !---- ORI ----!
2587
                   select case (nmol)
2588
                      case (-1)
2589
                         err_refcodes=.true.
2590
                         ERR_RefCodes_Mess="Option Not defined"
2591
                         return
2592

    
2593
                      case (0)
2594
                         do i=1,molcrys%n_mol
2595
                            do j=1,3
2596
                               if (molcrys%mol(i)%lOrient(j) /=0) then
2597
                                  nc=molcrys%mol(i)%lOrient(j)
2598
                                  call Delete_RefCodes(nc,molcrys)
2599
                               end if
2600
                            end do
2601
                         end do
2602

    
2603
                      case (1:)
2604
                         do j=1,3
2605
                            if (molcrys%mol(nmol)%lOrient(j) /=0) then
2606
                              nc=molcrys%mol(nmol)%lOrient(j)
2607
                               call Delete_RefCodes(nc,molcrys)
2608
                            end if
2609
                         end do
2610
                   end select
2611

    
2612
                case (8)
2613
                   !---- THE ----!
2614
                   !!! Not yet Implemented !!!!
2615
             end select
2616

    
2617
          !---- VARY Directive ----!
2618
          case ("var")
2619

    
2620
             select case (key)
2621
                case (0)
2622
                   !---- Nb must be different zero ----!
2623
                   select case (nb)
2624
                      case (0)
2625
                         err_refcodes=.true.
2626
                         ERR_RefCodes_Mess="Option Not defined"
2627
                         return
2628

    
2629
                      case ( 1:3)
2630
                         !---- X_, Y_, Z_ ----!
2631
                         if (na <= molcrys%n_free) then
2632
                            if (molcrys%atm(na)%lx(nb) ==0) then
2633
                               molcrys%atm(na)%mx(nb)=1.0
2634
                               call get_atompos_ctr(molcrys%atm(na)%x,   &
2635
                                                    molcrys%Spg,np_refi, &
2636
                                                    molcrys%atm(na)%lx,  &
2637
                                                    molcrys%atm(na)%mx)
2638
                               if (molcrys%atm(na)%lx(nb) == np_refi) then
2639
                                  V_Vec(np_refi)=molcrys%atm(na)%x(nb)
2640
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2641
                                  V_Bounds(1,np_refi)=xl
2642
                                  V_Bounds(2,np_refi)=xu
2643
                                  V_Bounds(3,np_refi)=xs
2644
                                  V_BCon(np_refi)=ic
2645
                                  V_list(np_refi)=na
2646
                               else
2647
                                  np_refi=np_refi-1
2648
                               end if
2649
                            end if
2650
                         else
2651
                            naa=na-molcrys%n_free
2652
                            do i=1,molcrys%n_mol
2653
                               if (naa > molcrys%mol(i)%natoms) then
2654
                                  naa=naa-molcrys%mol(i)%natoms
2655
                                  cycle
2656
                               end if
2657

    
2658
                               if (molcrys%mol(i)%lI_coor(nb,naa) ==0) then
2659
                                  molcrys%mol(i)%mI_coor(nb,naa)=1.0
2660
                                  call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa),  &
2661
                                                       molcrys%Spg, np_refi,          &
2662
                                                       molcrys%mol(i)%lI_coor(:,naa), &
2663
                                                       molcrys%mol(i)%mI_coor(:,naa))
2664
                                  if (molcrys%mol(i)%lI_coor(nb,naa) == np_refi) then
2665
                                     V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
2666
                                     V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2667
                                     V_Bounds(1,np_refi)=xl
2668
                                     V_Bounds(2,np_refi)=xu
2669
                                     V_Bounds(3,np_refi)=xs
2670
                                     V_BCon(np_refi)=ic
2671
                                     V_list(np_refi)=na
2672
                                  else
2673
                                     np_refi=np_refi-1
2674
                                  end if
2675
                               end if
2676
                            end do
2677
                         end if
2678

    
2679
                      case ( 4)
2680
                         !---- Biso_ ----!
2681
                         if (na <= molcrys%n_free) then
2682
                            if (molcrys%atm(na)%lbiso ==0) then
2683
                               np_refi=np_refi+1
2684
                               V_Vec(np_refi)=molcrys%atm(na)%biso
2685
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2686
                               molcrys%atm(na)%mbiso=1.0
2687
                               molcrys%atm(na)%lbiso=np_refi
2688
                               V_Bounds(1,np_refi)=xl
2689
                               V_Bounds(2,np_refi)=xu
2690
                               V_Bounds(3,np_refi)=xs
2691
                               V_BCon(np_refi)=ic
2692
                               V_list(np_refi)=na
2693
                            end if
2694
                         else
2695
                            naa=na-molcrys%n_free
2696
                            do i=1,molcrys%n_mol
2697
                               if (naa > molcrys%mol(i)%natoms) then
2698
                                  naa=naa-molcrys%mol(i)%natoms
2699
                                  cycle
2700
                               end if
2701
                               if (molcrys%mol(i)%lbiso(naa) ==0) then
2702
                                  np_refi=np_refi+1
2703
                                  V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
2704
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2705
                                  molcrys%mol(i)%mbiso(naa)=1.0
2706
                                  molcrys%mol(i)%lbiso(naa)=np_refi
2707
                                  V_Bounds(1,np_refi)=xl
2708
                                  V_Bounds(2,np_refi)=xu
2709
                                  V_Bounds(3,np_refi)=xs
2710
                                  V_BCon(np_refi)=ic
2711
                                  V_list(np_refi)=na
2712
                               end if
2713
                            end do
2714
                         end if
2715

    
2716
                      case ( 5)
2717
                         !---- Occ_ ----!
2718
                         if (na <= molcrys%n_free) then
2719
                            if (molcrys%atm(na)%locc ==0) then
2720
                               np_refi=np_refi+1
2721
                               V_Vec(np_refi)=molcrys%atm(na)%occ
2722
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2723
                               molcrys%atm(na)%mocc=1.0
2724
                               molcrys%atm(na)%locc=np_refi
2725
                               V_Bounds(1,np_refi)=xl
2726
                               V_Bounds(2,np_refi)=xu
2727
                               V_Bounds(3,np_refi)=xs
2728
                               V_BCon(np_refi)=ic
2729
                               V_list(np_refi)=na
2730
                            end if
2731
                         else
2732
                            naa=na-molcrys%n_free
2733
                            do i=1,molcrys%n_mol
2734
                               if (naa > molcrys%mol(i)%natoms) then
2735
                                  naa=naa-molcrys%mol(i)%natoms
2736
                                  cycle
2737
                               end if
2738
                               if (molcrys%mol(i)%locc(naa) ==0) then
2739
                                  np_refi=np_refi+1
2740
                                  V_Vec(np_refi)=molcrys%mol(i)%occ(naa)
2741
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2742
                                  molcrys%mol(i)%mocc(naa)=1.0
2743
                                  molcrys%mol(i)%locc(naa)=np_refi
2744
                                  V_Bounds(1,np_refi)=xl
2745
                                  V_Bounds(2,np_refi)=xu
2746
                                  V_Bounds(3,np_refi)=xs
2747
                                  V_BCon(np_refi)=ic
2748
                                  V_list(np_refi)=na
2749
                               end if
2750
                            end do
2751
                         end if
2752

    
2753
                      case ( 6:11)
2754
                         !---- B11_, ..., B23_ ----!
2755
                         if (na <= molcrys%n_free) then
2756
                            if (molcrys%atm(na)%lu(nb-5) ==0) then
2757
                               molcrys%atm(na)%mu(nb-5)=1.0
2758
                               call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
2759
                                                   np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
2760
                               if (molcrys%atm(na)%lu(nb-5) == np_refi) then
2761
                                  V_Vec(np_refi)=molcrys%atm(na)%u(nb-5)
2762
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2763
                                  V_Bounds(1,np_refi)=xl
2764
                                  V_Bounds(2,np_refi)=xu
2765
                                  V_Bounds(3,np_refi)=xs
2766
                                  V_BCon(np_refi)=ic
2767
                                  V_list(np_refi)=na
2768
                               else
2769
                                  np_refi=np_refi-1
2770
                               end if
2771
                            end if
2772
                         else
2773
                            err_refcodes=.true.
2774
                            ERR_RefCodes_Mess="Option Not defined"
2775
                            return
2776
                         end if
2777

    
2778
                      case (12)
2779
                         !---- Banis_ ----!
2780
                         if (na <= molcrys%n_free) then
2781
                            do j=1,6
2782
                               if (molcrys%atm(na)%lu(j) ==0) then
2783
                                  molcrys%atm(na)%mu(j)=1.0
2784
                                  call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
2785
                                                       np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
2786
                                  if (molcrys%atm(na)%lu(j) == np_refi) then
2787
                                     V_Vec(np_refi)=molcrys%atm(na)%u(j)
2788
                                     V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
2789
                                     V_Bounds(1,np_refi)=xl
2790
                                     V_Bounds(2,np_refi)=xu
2791
                                     V_Bounds(3,np_refi)=xs
2792
                                     V_BCon(np_refi)=ic
2793
                                     V_list(np_refi)=na
2794
                                  else
2795
                                     np_refi=np_refi-1
2796
                                  end if
2797
                               end if
2798
                            end do
2799
                         else
2800
                            err_refcodes=.true.
2801
                            ERR_RefCodes_Mess="Option Not defined"
2802
                            return
2803
                         end if
2804

    
2805
                      case (13:15)
2806
                         !---- Xc_, Yc_, Zc_ ----!
2807
                         select case (nmol)
2808
                            case (-1)
2809
                               err_refcodes=.true.
2810
                               ERR_RefCodes_Mess="Option Not defined"
2811
                               return
2812

    
2813
                            case (0)
2814
                               do i=1,molcrys%n_mol
2815
                                  write(unit=car,fmt="(i2)") i
2816
                                  car=adjustl(car)
2817
                                  if (molcrys%mol(i)%lxcentre(nb-12) ==0) then
2818
                                     np_refi=np_refi+1
2819
                                     V_Vec(np_refi)=molcrys%mol(i)%xcentre(nb-12)
2820
                                     V_Name(np_refi)=trim(code_nam(nb))//"Mol"//trim(car)
2821
                                     molcrys%mol(i)%mxcentre(nb-12)=1.0
2822
                                     molcrys%mol(i)%lxcentre(nb-12)=np_refi
2823
                                     V_Bounds(1,np_refi)=xl
2824
                                     V_Bounds(2,np_refi)=xu
2825
                                     V_Bounds(3,np_refi)=xs
2826
                                     V_BCon(np_refi)=ic
2827
                                     V_list(np_refi)=-i
2828
                                  end if
2829
                               end do
2830

    
2831
                            case (1:)
2832
                               write(unit=car,fmt="(i2)") nmol
2833
                               car=adjustl(car)
2834
                               if (molcrys%mol(nmol)%lxcentre(nb-12) ==0) then
2835
                                  np_refi=np_refi+1
2836
                                  V_Vec(np_refi)=molcrys%mol(nmol)%xcentre(nb-12)
2837
                                  V_Name(np_refi)=trim(code_nam(nb))//"entre_Mol"//trim(car)
2838
                                  molcrys%mol(nmol)%mxcentre(nb-12)=1.0
2839
                                  molcrys%mol(nmol)%lxcentre(nb-12)=np_refi
2840
                                  V_Bounds(1,np_refi)=xl
2841
                                  V_Bounds(2,np_refi)=xu
2842
                                  V_Bounds(3,np_refi)=xs
2843
                                  V_BCon(np_refi)=ic
2844
                                  V_list(np_refi)=-nmol
2845
                               end if
2846
                         end select
2847

    
2848
                      case (16:18)
2849
                         !---- Theta_, Phi_, Chi_ ----!
2850
                         select case (nmol)
2851
                            case (-1)
2852
                               err_refcodes=.true.
2853
                               ERR_RefCodes_Mess="Option Not defined"
2854
                               return
2855

    
2856
                            case (0)
2857
                               do i=1,molcrys%n_mol
2858
                                  write(unit=car,fmt="(i2)") i
2859
                                  car=adjustl(car)
2860
                                  if (molcrys%mol(i)%lOrient(nb-15) ==0) then
2861
                                     np_refi=np_refi+1
2862
                                     V_Vec(np_refi)=molcrys%mol(i)%Orient(nb-15)
2863
                                     V_Name(np_refi)=trim(code_nam(nb))//"Orient_Mol"//trim(car)
2864
                                     molcrys%mol(i)%mOrient(nb-15)=1.0
2865
                                     molcrys%mol(i)%lOrient(nb-15)=np_refi
2866
                                     V_Bounds(1,np_refi)=xl
2867
                                     V_Bounds(2,np_refi)=xu
2868
                                     V_Bounds(3,np_refi)=xs
2869
                                     V_BCon(np_refi)=ic
2870
                                     V_list(np_refi)=-i
2871
                                  end if
2872
                               end do
2873

    
2874
                            case (1:)
2875
                               write(unit=car,fmt="(i2)") nmol
2876
                               car=adjustl(car)
2877
                               if (molcrys%mol(nmol)%lOrient(nb-15) ==0) then
2878
                                  np_refi=np_refi+1
2879
                                  V_Vec(np_refi)=molcrys%mol(nmol)%Orient(nb-15)
2880
                                  V_Name(np_refi)=trim(code_nam(nb))//"Orient_Mol"//trim(car)
2881
                                  molcrys%mol(nmol)%mOrient(nb-15)=1.0
2882
                                  molcrys%mol(nmol)%lOrient(nb-15)=np_refi
2883
                                  V_Bounds(1,np_refi)=xl
2884
                                  V_Bounds(2,np_refi)=xu
2885
                                  V_Bounds(3,np_refi)=xs
2886
                                  V_BCon(np_refi)=ic
2887
                                  V_list(np_refi)=-nmol
2888
                               end if
2889
                         end select
2890

    
2891
                      case (19:21)
2892
                         !!! Not Yet Implemented !!!
2893

    
2894
                   end select ! nb
2895

    
2896
                case (1)
2897
                   !---- XYZ ----!
2898
                   if (na <= molcrys%n_free) then
2899
                      do j=1,3
2900
                         if (molcrys%atm(na)%lx(j) ==0) then
2901
                            molcrys%atm(na)%mx(j)=1.0
2902
                            call get_atompos_ctr(molcrys%atm(na)%x,   &
2903
                                                 molcrys%Spg,np_refi, &
2904
                                                 molcrys%atm(na)%lx,  &
2905
                                                 molcrys%atm(na)%mx)
2906
                            if (molcrys%atm(na)%lx(j) == np_refi) then
2907
                               V_Vec(np_refi)=molcrys%atm(na)%x(j)
2908
                               V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%atm(na)%lab)
2909
                               V_Bounds(1,np_refi)=xl
2910
                               V_Bounds(2,np_refi)=xu
2911
                               V_Bounds(3,np_refi)=xs
2912
                               V_BCon(np_refi)=ic
2913
                               V_list(np_refi)=na
2914
                            else
2915
                               np_refi=np_refi-1
2916
                            end if
2917
                         end if
2918
                      end do
2919
                   else
2920
                      naa=na-molcrys%n_free
2921
                      do i=1,molcrys%n_mol
2922
                         if (naa > molcrys%mol(i)%natoms) then
2923
                             naa=naa-molcrys%mol(i)%natoms
2924
                             cycle
2925
                         end if
2926
                         do j=1,3
2927
                            if (molcrys%mol(i)%lI_coor(j,naa) ==0) then
2928
                               molcrys%mol(i)%mI_coor(nb,naa)=1.0
2929
                               call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa), &
2930
                                                    molcrys%Spg,np_refi,   &
2931
                                                    molcrys%mol(i)%lI_coor(:,naa), &
2932
                                                    molcrys%mol(i)%mI_coor(:,naa))
2933
                               if (molcrys%mol(i)%lI_coor(j,naa) == np_refi) then
2934
                                  V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
2935
                                  V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%mol(i)%AtName(naa))
2936
                                  V_Bounds(1,np_refi)=xl
2937
                                  V_Bounds(2,np_refi)=xu
2938
                                  V_Bounds(3,np_refi)=xs
2939
                                  V_BCon(np_refi)=ic
2940
                                  V_list(np_refi)=na
2941
                               else
2942
                                  np_refi=np_refi-1
2943
                               end if
2944
                            end if
2945
                         end do
2946
                      end do
2947
                   end if
2948

    
2949
                case (2)
2950
                   !---- OCC ----!
2951
                   if (na <= molcrys%n_free) then
2952
                      if (molcrys%atm(na)%locc ==0) then
2953
                         np_refi=np_refi+1
2954
                         V_Vec(np_refi)=molcrys%atm(na)%occ
2955
                         V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%atm(na)%lab)
2956
                         molcrys%atm(na)%mocc=1.0
2957
                         molcrys%atm(na)%locc=np_refi
2958
                         V_Bounds(1,np_refi)=xl
2959
                         V_Bounds(2,np_refi)=xu
2960
                         V_Bounds(3,np_refi)=xs
2961
                         V_BCon(np_refi)=ic
2962
                         V_list(np_refi)=na
2963
                      end if
2964
                   else
2965
                      naa=na-molcrys%n_free
2966
                      do i=1,molcrys%n_mol
2967
                         if (naa > molcrys%mol(i)%natoms) then
2968
                            naa=naa-molcrys%mol(i)%natoms
2969
                            cycle
2970
                         end if
2971

    
2972
                         if (molcrys%mol(i)%locc(naa) ==0) then
2973
                            np_refi=np_refi+1
2974
                            V_Vec(np_refi)=molcrys%mol(i)%Occ(naa)
2975
                            V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%mol(i)%AtName(naa))
2976
                            molcrys%mol(i)%mocc(naa)=1.0
2977
                            molcrys%mol(i)%locc(naa)=np_refi
2978
                            V_Bounds(1,np_refi)=xl
2979
                            V_Bounds(2,np_refi)=xu
2980
                            V_Bounds(3,np_refi)=xs
2981
                            V_BCon(np_refi)=ic
2982
                            V_list(np_refi)=na
2983
                         end if
2984
                      end do
2985
                   end if
2986

    
2987
                case (3)
2988
                   !---- BIS ----!
2989
                   if (na <= molcrys%n_free) then
2990
                      if (molcrys%atm(na)%lbiso ==0) then
2991
                         np_refi=np_refi+1
2992
                         V_Vec(np_refi)=molcrys%atm(na)%biso
2993
                         V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%atm(na)%lab)
2994
                         molcrys%atm(na)%mbiso=1.0
2995
                         molcrys%atm(na)%lbiso=np_refi
2996
                         V_Bounds(1,np_refi)=xl
2997
                         V_Bounds(2,np_refi)=xu
2998
                         V_Bounds(3,np_refi)=xs
2999
                         V_BCon(np_refi)=ic
3000
                         V_list(np_refi)=na
3001
                      end if
3002
                   else
3003
                      naa=na-molcrys%n_free
3004
                      do i=1,molcrys%n_mol
3005
                         if (naa > molcrys%mol(i)%natoms) then
3006
                            naa=naa-molcrys%mol(i)%natoms
3007
                            cycle
3008
                         end if
3009
                         if (molcrys%mol(i)%lbiso(naa) ==0) then
3010
                            np_refi=np_refi+1
3011
                            V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
3012
                            V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%mol(i)%AtName(naa))
3013
                            molcrys%mol(i)%mbiso(naa)=1.0
3014
                            molcrys%mol(i)%lbiso(naa)=np_refi
3015
                            V_Bounds(1,np_refi)=xl
3016
                            V_Bounds(2,np_refi)=xu
3017
                            V_Bounds(3,np_refi)=xs
3018
                            V_BCon(np_refi)=ic
3019
                            V_list(np_refi)=na
3020
                         end if
3021
                      end do
3022
                   end if
3023

    
3024
                case (4)
3025
                   !---- BAN ----!
3026
                   if (na <= molcrys%n_free) then
3027
                      do j=1,6
3028
                         if (molcrys%atm(na)%lu(j) ==0) then
3029
                            molcrys%atm(na)%mu(j)=1.0
3030
                            call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
3031
                                                 np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
3032
                            if (molcrys%atm(na)%lu(j) == np_refi) then
3033
                               V_Vec(np_refi)=molcrys%atm(na)%u(j)
3034
                               V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
3035
                               V_Bounds(1,np_refi)=xl
3036
                               V_Bounds(2,np_refi)=xu
3037
                               V_Bounds(3,np_refi)=xs
3038
                               V_BCon(np_refi)=ic
3039
                               V_list(np_refi)=na
3040
                            else
3041
                               np_refi=np_refi-1
3042
                            end if
3043
                         end if
3044
                      end do
3045
                   else
3046
                      err_refcodes=.true.
3047
                      ERR_RefCodes_Mess="Option Not defined"
3048
                      return
3049
                   end if
3050

    
3051
                case (5)
3052
                   !---- ALL ----!
3053
                   if (na <= molcrys%n_free) then
3054
                      do j=1,3
3055
                         if (molcrys%atm(na)%lx(j) ==0) then
3056
                            molcrys%atm(na)%mx(j)=1.0
3057
                            call get_atompos_ctr(molcrys%atm(na)%x,   &
3058
                                                 molcrys%Spg,np_refi, &
3059
                                                 molcrys%atm(na)%lx,  &
3060
                                                 molcrys%atm(na)%mx)
3061
                            if (molcrys%atm(na)%lx(j) == np_refi) then
3062
                               V_Vec(np_refi)=molcrys%atm(na)%x(j)
3063
                               V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%atm(na)%lab)
3064
                               V_Bounds(1,np_refi)=xl
3065
                               V_Bounds(2,np_refi)=xu
3066
                               V_Bounds(3,np_refi)=xs
3067
                               V_BCon(np_refi)=ic
3068
                               V_list(np_refi)=na
3069
                            else
3070
                               np_refi=np_refi-1
3071
                            end if
3072
                         end if
3073
                      end do
3074
                   else
3075
                      naa=na-molcrys%n_free
3076
                      do i=1,molcrys%n_mol
3077
                         if (naa > molcrys%mol(i)%natoms) then
3078
                             naa=naa-molcrys%mol(i)%natoms
3079
                             cycle
3080
                         end if
3081

    
3082
                         do j=1,3
3083
                            if (molcrys%mol(i)%lI_coor(j,naa) ==0) then
3084
                               molcrys%mol(i)%mI_coor(nb,naa)=1.0
3085
                               call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa), &
3086
                                                    molcrys%Spg,np_refi,   &
3087
                                                    molcrys%mol(i)%lI_coor(:,naa), &
3088
                                                    molcrys%mol(i)%mI_coor(:,naa))
3089
                               if (molcrys%mol(i)%lI_coor(j,naa) == np_refi) then
3090
                                  V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
3091
                                  V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%mol(i)%AtName(naa))
3092
                                  V_Bounds(1,np_refi)=xl
3093
                                  V_Bounds(2,np_refi)=xu
3094
                                  V_Bounds(3,np_refi)=xs
3095
                                  V_BCon(np_refi)=ic
3096
                                  V_list(np_refi)=na
3097
                               else
3098
                                  np_refi=np_refi-1
3099
                               end if
3100
                            end if
3101
                         end do
3102
                      end do
3103
                   end if
3104

    
3105
                   if (na <= molcrys%n_free) then
3106
                      if (molcrys%atm(na)%locc ==0) then
3107
                         np_refi=np_refi+1
3108
                         V_Vec(np_refi)=molcrys%atm(na)%occ
3109
                         V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%atm(na)%lab)
3110
                         molcrys%atm(na)%mocc=1.0
3111
                         molcrys%atm(na)%locc=np_refi
3112
                         V_Bounds(1,np_refi)=xl
3113
                         V_Bounds(2,np_refi)=xu
3114
                         V_Bounds(3,np_refi)=xs
3115
                         V_BCon(np_refi)=ic
3116
                         V_list(np_refi)=na
3117
                      end if
3118
                   else
3119
                      naa=na-molcrys%n_free
3120
                      do i=1,molcrys%n_mol
3121
                         if (naa > molcrys%mol(i)%natoms) then
3122
                            naa=naa-molcrys%mol(i)%natoms
3123
                            cycle
3124
                         end if
3125
                         if (molcrys%mol(i)%locc(naa) ==0) then
3126
                            np_refi=np_refi+1
3127
                            V_Vec(np_refi)=molcrys%mol(i)%Occ(naa)
3128
                            V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%mol(i)%AtName(naa))
3129
                            molcrys%mol(i)%mocc(naa)=1.0
3130
                            molcrys%mol(i)%locc(naa)=np_refi
3131
                            V_Bounds(1,np_refi)=xl
3132
                            V_Bounds(2,np_refi)=xu
3133
                            V_Bounds(3,np_refi)=xs
3134
                            V_BCon(np_refi)=ic
3135
                            V_list(np_refi)=na
3136
                         end if
3137
                      end do
3138
                   end if
3139

    
3140
                   if (na <= molcrys%n_free) then
3141
                      if (molcrys%atm(na)%lbiso ==0) then
3142
                         np_refi=np_refi+1
3143
                         V_Vec(np_refi)=molcrys%atm(na)%biso
3144
                         V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%atm(na)%lab)
3145
                         molcrys%atm(na)%mbiso=1.0
3146
                         molcrys%atm(na)%lbiso=np_refi
3147
                         V_Bounds(1,np_refi)=xl
3148
                         V_Bounds(2,np_refi)=xu
3149
                         V_Bounds(3,np_refi)=xs
3150
                         V_BCon(np_refi)=ic
3151
                         V_list(np_refi)=na
3152
                      end if
3153
                   else
3154
                      naa=na-molcrys%n_free
3155
                      do i=1,molcrys%n_mol
3156
                         if (naa > molcrys%mol(i)%natoms) then
3157
                            naa=naa-molcrys%mol(i)%natoms
3158
                            cycle
3159
                         end if
3160

    
3161
                         if (molcrys%mol(i)%lbiso(naa) ==0) then
3162
                            np_refi=np_refi+1
3163
                            V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
3164
                            V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%mol(i)%AtName(naa))
3165
                            molcrys%mol(i)%mbiso(naa)=1.0
3166
                            molcrys%mol(i)%lbiso(naa)=np_refi
3167
                            V_Bounds(1,np_refi)=xl
3168
                            V_Bounds(2,np_refi)=xu
3169
                            V_Bounds(3,np_refi)=xs
3170
                            V_BCon(np_refi)=ic
3171
                            V_list(np_refi)=na
3172
                         end if
3173
                      end do
3174
                   end if
3175

    
3176
                   if (na <= molcrys%n_free) then
3177
                      do j=1,6
3178
                         if (molcrys%atm(na)%lu(j) ==0) then
3179
                            molcrys%atm(na)%mu(j)=1.0
3180
                            call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
3181
                                                 np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
3182
                            if (molcrys%atm(na)%lu(j) == np_refi) then
3183
                               V_Vec(np_refi)=molcrys%atm(na)%u(j)
3184
                               V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
3185
                               V_Bounds(1,np_refi)=xl
3186
                               V_Bounds(2,np_refi)=xu
3187
                               V_Bounds(3,np_refi)=xs
3188
                               V_BCon(np_refi)=ic
3189
                               V_list(np_refi)=na
3190
                            else
3191
                               np_refi=np_refi-1
3192
                            end if
3193
                         end if
3194
                      end do
3195
                   end if
3196

    
3197
                   select case (nmol)
3198
                      case (-1)
3199
                         err_refcodes=.true.
3200
                         ERR_RefCodes_Mess="Option Not defined"
3201
                         return
3202

    
3203
                      case (0)
3204
                         do i=1,molcrys%n_mol
3205
                            write(unit=car,fmt="(i2)") i
3206
                            car=adjustl(car)
3207
                            do j=1,3
3208
                               if (molcrys%mol(i)%lxcentre(j) ==0) then
3209
                                  np_refi=np_refi+1
3210
                                  V_Vec(np_refi)=molcrys%mol(i)%xcentre(j)
3211
                                  V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3212
                                  molcrys%mol(i)%mxcentre(j)=1.0
3213
                                  molcrys%mol(i)%lxcentre(j)=np_refi
3214
                                  V_Bounds(1,np_refi)=xl
3215
                                  V_Bounds(2,np_refi)=xu
3216
                                  V_Bounds(3,np_refi)=xs
3217
                                  V_BCon(np_refi)=ic
3218
                                  V_list(np_refi)=-i
3219
                               end if
3220
                            end do
3221
                         end do
3222

    
3223
                      case (1:)
3224
                         write(unit=car,fmt="(i2)") nmol
3225
                         car=adjustl(car)
3226
                         do j=1,3
3227
                            if (molcrys%mol(nmol)%lxcentre(j) ==0) then
3228
                               np_refi=np_refi+1
3229
                               V_Vec(np_refi)=molcrys%mol(nmol)%xcentre(j)
3230
                               V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3231
                               molcrys%mol(nmol)%mxcentre(j)=1.0
3232
                               molcrys%mol(nmol)%lxcentre(j)=np_refi
3233
                               V_Bounds(1,np_refi)=xl
3234
                               V_Bounds(2,np_refi)=xu
3235
                               V_Bounds(3,np_refi)=xs
3236
                               V_BCon(np_refi)=ic
3237
                               V_list(np_refi)=-nmol
3238
                            end if
3239
                         end do
3240
                   end select
3241

    
3242
                   select case (nmol)
3243
                      case (-1)
3244
                         err_refcodes=.true.
3245
                         ERR_RefCodes_Mess="Option Not defined"
3246
                         return
3247

    
3248
                      case (0)
3249
                         do i=1,molcrys%n_mol
3250
                            write(unit=car,fmt="(i2)") i
3251
                            car=adjustl(car)
3252
                            do j=1,3
3253
                               if (molcrys%mol(i)%lOrient(j) ==0) then
3254
                                  np_refi=np_refi+1
3255
                                  V_Vec(np_refi)=molcrys%mol(i)%Orient(j)
3256
                                  V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3257
                                  molcrys%mol(i)%mOrient(j)=1.0
3258
                                  molcrys%mol(i)%lOrient(j)=np_refi
3259
                                  V_Bounds(1,np_refi)=xl
3260
                                  V_Bounds(2,np_refi)=xu
3261
                                  V_Bounds(3,np_refi)=xs
3262
                                  V_BCon(np_refi)=ic
3263
                                  V_list(np_refi)=-i
3264
                               end if
3265
                            end do
3266
                         end do
3267

    
3268
                      case (1:)
3269
                         write(unit=car,fmt="(i2)") nmol
3270
                         car=adjustl(car)
3271
                         do j=1,3
3272
                            if (molcrys%mol(nmol)%lOrient(j) ==0) then
3273
                               np_refi=np_refi+1
3274
                               V_Vec(np_refi)=molcrys%mol(nmol)%Orient(j)
3275
                               V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3276
                               molcrys%mol(nmol)%mOrient(j)=1.0
3277
                               molcrys%mol(nmol)%lOrient(j)=np_refi
3278
                               V_Bounds(1,np_refi)=xl
3279
                               V_Bounds(2,np_refi)=xu
3280
                               V_Bounds(3,np_refi)=xs
3281
                               V_BCon(np_refi)=ic
3282
                               V_list(np_refi)=-nmol
3283
                            end if
3284
                         end do
3285
                   end select
3286

    
3287
                case (6)
3288
                   !---- CEN ----!
3289
                   select case (nmol)
3290
                      case (-1)
3291
                         err_refcodes=.true.
3292
                         ERR_RefCodes_Mess="Option Not defined"
3293
                         return
3294

    
3295
                      case (0)
3296
                         do i=1,molcrys%n_mol
3297
                            write(unit=car,fmt="(i2)") i
3298
                            car=adjustl(car)
3299
                            do j=1,3
3300
                               if (molcrys%mol(i)%lxcentre(j) ==0) then
3301
                                  np_refi=np_refi+1
3302
                                  V_Vec(np_refi)=molcrys%mol(i)%xcentre(j)
3303
                                  V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3304
                                  molcrys%mol(i)%mxcentre(j)=1.0
3305
                                  molcrys%mol(i)%lxcentre(j)=np_refi
3306
                                  V_Bounds(1,np_refi)=xl
3307
                                  V_Bounds(2,np_refi)=xu
3308
                                  V_Bounds(3,np_refi)=xs
3309
                                  V_BCon(np_refi)=ic
3310
                                  V_list(np_refi)=-i
3311
                               end if
3312
                            end do
3313
                         end do
3314

    
3315
                      case (1:)
3316
                         write(unit=car,fmt="(i2)") nmol
3317
                         car=adjustl(car)
3318
                         do j=1,3
3319
                            if (molcrys%mol(nmol)%lxcentre(j) ==0) then
3320
                               np_refi=np_refi+1
3321
                               V_Vec(np_refi)=molcrys%mol(nmol)%xcentre(j)
3322
                               V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3323
                               molcrys%mol(nmol)%mxcentre(j)=1.0
3324
                               molcrys%mol(nmol)%lxcentre(j)=np_refi
3325
                               V_Bounds(1,np_refi)=xl
3326
                               V_Bounds(2,np_refi)=xu
3327
                               V_Bounds(3,np_refi)=xs
3328
                               V_BCon(np_refi)=ic
3329
                               V_list(np_refi)=-nmol
3330
                            end if
3331
                         end do
3332
                   end select
3333

    
3334
                case (7)
3335
                   !---- ORI  ----!
3336
                   select case (nmol)
3337
                      case (-1)
3338
                         err_refcodes=.true.
3339
                         ERR_RefCodes_Mess="Option Not defined"
3340
                         return
3341

    
3342
                      case (0)
3343
                         do i=1,molcrys%n_mol
3344
                            write(unit=car,fmt="(i2)") i
3345
                            car=adjustl(car)
3346
                            do j=1,3
3347
                               if (molcrys%mol(i)%lOrient(j) ==0) then
3348
                                  np_refi=np_refi+1
3349
                                  V_Vec(np_refi)=molcrys%mol(i)%Orient(j)
3350
                                  V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3351
                                  molcrys%mol(i)%mOrient(j)=1.0
3352
                                  molcrys%mol(i)%lOrient(j)=np_refi
3353
                                  V_Bounds(1,np_refi)=xl
3354
                                  V_Bounds(2,np_refi)=xu
3355
                                  V_Bounds(3,np_refi)=xs
3356
                                  V_BCon(np_refi)=ic
3357
                                  V_list(np_refi)=-i
3358
                               end if
3359
                            end do
3360
                         end do
3361

    
3362
                      case (1:)
3363
                         write(unit=car,fmt="(i2)") nmol
3364
                         car=adjustl(car)
3365
                         do j=1,3
3366
                            if (molcrys%mol(nmol)%lOrient(j) ==0) then
3367
                               np_refi=np_refi+1
3368
                               V_Vec(np_refi)=molcrys%mol(nmol)%Orient(j)
3369
                               V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3370
                               molcrys%mol(nmol)%mOrient(j)=1.0
3371
                               molcrys%mol(nmol)%lOrient(j)=np_refi
3372
                               V_Bounds(1,np_refi)=xl
3373
                               V_Bounds(2,np_refi)=xu
3374
                               V_Bounds(3,np_refi)=xs
3375
                               V_BCon(np_refi)=ic
3376
                               V_list(np_refi)=-nmol
3377
                            end if
3378
                         end do
3379
                   end select
3380

    
3381
                case (8)
3382
                   !!! Not yet implemented !!!
3383

    
3384
             end select
3385
       end select
3386

    
3387
       return
3388
    End Subroutine Fill_RefCodes_Molcrys
3389

    
3390
    !!--++
3391
    !!--++ Subroutine Fill_RefCodes_Molec(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molec,Spg)
3392
    !!--++    integer,                      intent(in)     :: Key
3393
    !!--++    character(len=*),             intent(in)     :: Dire
3394
    !!--++    integer,                      intent(in)     :: Na
3395
    !!--++    integer,                      intent(in)     :: Nb
3396
    !!--++    real(kind=cp),                intent(in)     :: Xl
3397
    !!--++    real(kind=cp),                intent(in)     :: Xu
3398
    !!--++    real(kind=cp),                intent(in)     :: Xs
3399
    !!--++    integer,                      intent(in)     :: Ic
3400
    !!--++    type(molecule_type),          intent(in out) :: Molec
3401
    !!--++    type(space_group_type),       intent(in)     :: Spg
3402
    !!--++
3403
    !!--++ Overloaded
3404
    !!--++ Write on Vectors the Information for Molecule_Type
3405
    !!--++
3406
    !!--++ Update: March - 2005
3407
    !!
3408
    Subroutine Fill_RefCodes_Molec(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molec,Spg)
3409
       !---- Arguments ----!
3410
       integer,                      intent(in)     :: Key
3411
       character(len=*),             intent(in)     :: Dire
3412
       integer,                      intent(in)     :: Na
3413
       integer,                      intent(in)     :: Nb
3414
       real(kind=cp),                intent(in)     :: Xl
3415
       real(kind=cp),                intent(in)     :: Xu
3416
       real(kind=cp),                intent(in)     :: Xs
3417
       integer,                      intent(in)     :: Ic
3418
       type(molecule_type),          intent(in out) :: Molec
3419
       type(space_group_type),       intent(in)     :: Spg
3420

    
3421
       !---- Local variables ----!
3422
       integer :: j, nc
3423

    
3424
       call init_err_refcodes()
3425
       if (Na <= 0) then
3426
          err_refcodes=.true.
3427
          ERR_RefCodes_Mess="Number of atom no defined"
3428
          return
3429
       end if
3430

    
3431
       select case (dire)
3432
          !---- FIX Directive ----!
3433
          case ("fix")
3434

    
3435
             select case (key)
3436
                case (0)
3437

    
3438
                   !---- nb must be different zero ----!
3439
                   select case (nb)
3440
                      case (0)
3441
                         err_refcodes=.true.
3442
                         ERR_RefCodes_Mess="Option not defined"
3443
                         return
3444

    
3445
                      case ( 1:3)
3446
                         !---- X_, Y_, Z_----!
3447
                         if (molec%lI_coor(nb,na) /=0) then
3448
                            nc=molec%lI_coor(nb,na)
3449
                            call Delete_RefCodes(nc,molec)
3450
                         end if
3451

    
3452
                      case ( 4)
3453
                         !---- Biso_ ----!
3454
                         if (molec%lbiso(na) /=0) then
3455
                            nc=molec%lbiso(na)
3456
                            call Delete_RefCodes(nc,molec)
3457
                         end if
3458

    
3459
                      case ( 5)
3460
                         !---- Occ_ ----!
3461
                         if (molec%locc(na) /=0) then
3462
                            nc=molec%locc(na)
3463
                            call Delete_RefCodes(nc,molec)
3464
                         end if
3465

    
3466
                      case ( 6:12)
3467
                         !---- B11_, ..., B23_ ----!
3468
                         err_refcodes=.true.
3469
                         ERR_RefCodes_Mess="Option not defined"
3470
                         return
3471

    
3472
                      case (13:15)
3473
                         !---- Xc_, Yc_, Zc_ ----!
3474
                         if (molec%lxcentre(nb-12) /=0) then
3475
                            nc=molec%lxcentre(nb-12)
3476
                            call Delete_RefCodes(nc,molec)
3477
                         end if
3478

    
3479
                      case (16:18)
3480
                         !---- Theta_, Phi_, Chi_ ----!
3481
                         if (molec%lOrient(nb-15) /=0) then
3482
                            nc=molec%lOrient(nb-15)
3483
                            call Delete_RefCodes(nc,molec)
3484
                         end if
3485

    
3486
                      case (19:21)
3487
                         !!! Not yet implement !!!
3488

    
3489
                   end select ! nb
3490

    
3491
                case (1)
3492
                   !---- XYZ ----!
3493
                   do j=1,3
3494
                      if (molec%lI_coor(j,na) /=0) then
3495
                         nc=molec%lI_coor(j,na)
3496
                         call Delete_RefCodes(nc,molec)
3497
                      end if
3498
                   end do
3499

    
3500
                case (2)
3501
                   !---- OCC ----!
3502
                   if (molec%locc(na) /=0) then
3503
                      nc=molec%locc(na)
3504
                      call Delete_RefCodes(nc,molec)
3505
                   end if
3506

    
3507
                case (3)
3508
                   !---- BIS ----!
3509
                   if (molec%lbiso(na) /=0) then
3510
                      nc=molec%lbiso(na)
3511
                      call Delete_RefCodes(nc,molec)
3512
                   end if
3513

    
3514
                case (4)
3515
                   !---- BAN ----!
3516
                   err_refcodes=.true.
3517
                   ERR_RefCodes_Mess="Option not defined"
3518
                   return
3519

    
3520
                case (5)
3521
                  !---- ALL ----!
3522
                  do j=1,3
3523
                     if (molec%lI_coor(j,na) /=0) then
3524
                        nc=molec%lI_coor(j,na)
3525
                        call Delete_RefCodes(nc,molec)
3526
                     end if
3527
                  end do
3528
                  if (molec%locc(na) /=0) then
3529
                     nc=molec%locc(na)
3530
                     call Delete_RefCodes(nc,molec)
3531
                  end if
3532
                  if (molec%lbiso(na) /=0) then
3533
                     nc=molec%lbiso(na)
3534
                     call Delete_RefCodes(nc,molec)
3535
                  end if
3536
                  do j=1,3
3537
                     if (molec%lxcentre(j) /=0) then
3538
                        nc=molec%lxcentre(j)
3539
                        call Delete_RefCodes(nc,molec)
3540
                     end if
3541
                  end do
3542
                  do j=1,3
3543
                     if (molec%lOrient(j) /=0) then
3544
                        nc=molec%lOrient(j)
3545
                        call Delete_RefCodes(nc,molec)
3546
                     end if
3547
                  end do
3548
                  !!! Falta THermal TLS !!!
3549

    
3550
               case (6)
3551
                  !---- CEN ----!
3552
                  do j=1,3
3553
                     if (molec%lxcentre(j) /=0) then
3554
                        nc=molec%lxcentre(j)
3555
                        call Delete_RefCodes(nc,molec)
3556
                     end if
3557
                  end do
3558

    
3559
               case (7)
3560
                  !---- ORI ----!
3561
                  do j=1,3
3562
                     if (molec%lOrient(j) /=0) then
3563
                        nc=molec%lOrient(j)
3564
                        call Delete_RefCodes(nc,molec)
3565
                     end if
3566
                  end do
3567

    
3568
               case (8)
3569
                  !---- THE ----!
3570
                  !!! Not Yet Implemented !!!
3571
             end select
3572

    
3573
          !---- VARY Directive ----!
3574
          case ("var")
3575

    
3576
             select case (key)
3577
                case (0)
3578

    
3579
                   !---- Nb must be different zero ----!
3580
                   select case (nb)
3581
                      case (0)
3582
                         err_refcodes=.true.
3583
                         ERR_RefCodes_Mess="Option not defined"
3584
                         return
3585

    
3586
                      case ( 1:3)
3587
                         !--- X_, Y_, Z_ ----!
3588
                         if (molec%lI_coor(nb,na) ==0) then
3589
                            molec%mI_coor(nb,na)=1.0
3590
                            call get_atompos_ctr(molec%I_Coor(:,na),  &
3591
                                                 Spg, np_refi,  &
3592
                                                 molec%lI_coor(:,na), &
3593
                                                 molec%mI_coor(:,na))
3594
                            if (molec%lI_coor(nb,na) == np_refi) then
3595
                               V_Vec(np_refi)=molec%I_Coor(nb,na)
3596
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3597
                               V_Bounds(1,np_refi)=xl
3598
                               V_Bounds(2,np_refi)=xu
3599
                               V_Bounds(3,np_refi)=xs
3600
                               V_BCon(np_refi)=ic
3601
                               V_list(np_refi)=na
3602
                            else
3603
                               np_refi=np_refi-1
3604
                            end if
3605
                         end if
3606

    
3607
                      case ( 4)
3608
                         !---- Biso_ ----!
3609
                         if (molec%lbiso(na) ==0) then
3610
                            np_refi=np_refi+1
3611
                            V_Vec(np_refi)=molec%biso(na)
3612
                            V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3613
                            molec%mbiso(na)=1.0
3614
                            molec%lbiso(na)=np_refi
3615
                            V_Bounds(1,np_refi)=xl
3616
                            V_Bounds(2,np_refi)=xu
3617
                            V_Bounds(3,np_refi)=xs
3618
                            V_BCon(np_refi)=ic
3619
                            V_list(np_refi)=na
3620
                         end if
3621

    
3622
                      case ( 5)
3623
                         !---- Occ_ ----!
3624
                         if (molec%locc(na) ==0) then
3625
                            np_refi=np_refi+1
3626
                            V_Vec(np_refi)=molec%occ(na)
3627
                            V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3628
                            molec%mocc(na)=1.0
3629
                            molec%locc(na)=np_refi
3630
                            V_Bounds(1,np_refi)=xl
3631
                            V_Bounds(2,np_refi)=xu
3632
                            V_Bounds(3,np_refi)=xs
3633
                            V_BCon(np_refi)=ic
3634
                            V_list(np_refi)=na
3635
                         end if
3636

    
3637
                      case ( 6:12)
3638
                         !---- B11_, ..., B23_ ----!
3639
                         err_refcodes=.true.
3640
                         ERR_RefCodes_Mess="Option not defined"
3641
                         return
3642

    
3643
                      case (13:15)
3644
                         !---- Xc_, Yc_, Zc_ ----!
3645
                         if (molec%lxcentre(nb-12) ==0) then
3646
                            np_refi=np_refi+1
3647
                            V_Vec(np_refi)=molec%xcentre(nb-12)
3648
                            V_Name(np_refi)=trim(code_nam(nb))//"Mol"
3649
                            molec%mxcentre(nb-12)=1.0
3650
                            molec%lxcentre(nb-12)=np_refi
3651
                            V_Bounds(1,np_refi)=xl
3652
                            V_Bounds(2,np_refi)=xu
3653
                            V_Bounds(3,np_refi)=xs
3654
                            V_BCon(np_refi)=ic
3655
                            V_list(np_refi)=0
3656
                         end if
3657

    
3658
                      case (16:18)
3659
                         !---- Theta_, Phi_, Chi_ ----!
3660
                         if (molec%lOrient(nb-15) ==0) then
3661
                            np_refi=np_refi+1
3662
                            V_Vec(np_refi)=molec%orient(nb-15)
3663
                            V_Name(np_refi)=trim(code_nam(nb))//"Mol"
3664
                            molec%mOrient(nb-15)=1.0
3665
                            molec%lOrient(nb-15)=np_refi
3666
                            V_Bounds(1,np_refi)=xl
3667
                            V_Bounds(2,np_refi)=xu
3668
                            V_Bounds(3,np_refi)=xs
3669
                            V_BCon(np_refi)=0
3670
                            V_list(np_refi)=0
3671
                         end if
3672

    
3673
                      case (19:21)
3674
                         !!! Not yet implement !!!
3675

    
3676
                   end select ! nb
3677

    
3678
                case (1)
3679
                   !---- XYZ ----!
3680
                   do j=1,3
3681
                      if (molec%lI_coor(j,na) ==0) then
3682
                         molec%mI_coor(j,na)=1.0
3683
                         call get_atompos_ctr(molec%I_Coor(:,na),  &
3684
                                              Spg, np_refi,  &
3685
                                              molec%lI_coor(:,na), &
3686
                                              molec%mI_coor(:,na))
3687
                         if (molec%lI_coor(j,na) == np_refi) then
3688
                            V_Vec(np_refi)=molec%I_Coor(j,na)
3689
                            V_Name(np_refi)=trim(code_nam(j))//trim(molec%AtName(na))
3690
                            V_Bounds(1,np_refi)=xl
3691
                            V_Bounds(2,np_refi)=xu
3692
                            V_Bounds(3,np_refi)=xs
3693
                            V_BCon(np_refi)=ic
3694
                            V_list(np_refi)=na
3695
                         else
3696
                            np_refi=np_refi-1
3697
                         end if
3698
                      end if
3699
                   end do
3700

    
3701
                case (2)
3702
                   !---- OCC ----!
3703
                   if (molec%locc(na) ==0) then
3704
                      np_refi=np_refi+1
3705
                      V_Vec(np_refi)=molec%occ(na)
3706
                      V_Name(np_refi)=trim(code_nam(5))//trim(molec%AtName(na))
3707
                      molec%mocc(na)=1.0
3708
                      molec%locc(na)=np_refi
3709
                      V_Bounds(1,np_refi)=xl
3710
                      V_Bounds(2,np_refi)=xu
3711
                      V_Bounds(3,np_refi)=xs
3712
                      V_BCon(np_refi)=ic
3713
                      V_list(np_refi)=na
3714
                   end if
3715

    
3716
                case (3)
3717
                   !---- BIS ----!
3718
                   if (molec%lbiso(na) ==0) then
3719
                      np_refi=np_refi+1
3720
                      V_Vec(np_refi)=molec%biso(na)
3721
                      V_Name(np_refi)=trim(code_nam(4))//trim(molec%AtName(na))
3722
                      molec%mbiso(na)=1.0
3723
                      molec%lbiso(na)=np_refi
3724
                      V_Bounds(1,np_refi)=xl
3725
                      V_Bounds(2,np_refi)=xu
3726
                      V_Bounds(3,np_refi)=xs
3727
                      V_BCon(np_refi)=ic
3728
                      V_list(np_refi)=na
3729
                   end if
3730

    
3731
                case (4)
3732
                   !---- BAN ----!
3733
                   err_refcodes=.true.
3734
                   ERR_RefCodes_Mess="Option not defined"
3735
                   return
3736

    
3737
                case (5)
3738
                   !---- ALL ----!
3739
                   do j=1,3
3740
                      if (molec%lI_coor(j,na) ==0) then
3741
                         molec%mI_coor(j,na)=1.0
3742
                         call get_atompos_ctr(molec%I_Coor(:,na),  &
3743
                                              Spg, np_refi,  &
3744
                                              molec%lI_coor(:,na), &
3745
                                              molec%mI_coor(:,na))
3746
                         if (molec%lI_coor(j,na) == np_refi) then
3747
                            V_Vec(np_refi)=molec%I_Coor(j,na)
3748
                            V_Name(np_refi)=trim(code_nam(j))//trim(molec%AtName(na))
3749
                            V_Bounds(1,np_refi)=xl
3750
                            V_Bounds(2,np_refi)=xu
3751
                            V_Bounds(3,np_refi)=xs
3752
                            V_BCon(np_refi)=ic
3753
                            V_list(np_refi)=na
3754
                         else
3755
                            np_refi=np_refi-1
3756
                         end if
3757
                      end if
3758
                   end do
3759
                   if (molec%locc(na) ==0) then
3760
                      np_refi=np_refi+1
3761
                      V_Vec(np_refi)=molec%occ(na)
3762
                      V_Name(np_refi)=trim(code_nam(5))//trim(molec%AtName(na))
3763
                      molec%mocc(na)=1.0
3764
                      molec%locc(na)=np_refi
3765
                      V_Bounds(1,np_refi)=xl
3766
                      V_Bounds(2,np_refi)=xu
3767
                      V_Bounds(3,np_refi)=xs
3768
                      V_BCon(np_refi)=ic
3769
                      V_list(np_refi)=na
3770
                   end if
3771
                   if (molec%lbiso(na) ==0) then
3772
                      np_refi=np_refi+1
3773
                      V_Vec(np_refi)=molec%biso(na)
3774
                      V_Name(np_refi)=trim(code_nam(4))//trim(molec%AtName(na))
3775
                      molec%mbiso(na)=1.0
3776
                      molec%lbiso(na)=np_refi
3777
                      V_Bounds(1,np_refi)=xl
3778
                      V_Bounds(2,np_refi)=xu
3779
                      V_Bounds(3,np_refi)=xs
3780
                      V_BCon(np_refi)=ic
3781
                      V_list(np_refi)=na
3782
                   end if
3783
                   do j=1,3
3784
                      if (molec%lxcentre(j) ==0) then
3785
                         np_refi=np_refi+1
3786
                         V_Vec(np_refi)=molec%xcentre(j)
3787
                         V_name(np_refi)=trim(code_nam(12+j))//"entre_Mol"
3788
                         molec%mxcentre(j)=1.0
3789
                         molec%lxcentre(j)=np_refi
3790
                         V_Bounds(1,np_refi)=xl
3791
                         V_Bounds(2,np_refi)=xu
3792
                         V_Bounds(3,np_refi)=xs
3793
                         V_BCon(np_refi)=ic
3794
                         V_list(np_refi)=0
3795
                      end if
3796
                   end do
3797
                   do j=1,3
3798
                      if (molec%lOrient(j) ==0) then
3799
                         np_refi=np_refi+1
3800
                         V_Vec(np_refi)=molec%orient(j)
3801
                         V_name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"
3802
                         molec%mOrient(j)=1.0
3803
                         molec%lOrient(j)=np_refi
3804
                         V_Bounds(1,np_refi)=xl
3805
                         V_Bounds(2,np_refi)=xu
3806
                         V_Bounds(3,np_refi)=xs
3807
                         V_BCon(np_refi)=0
3808
                         V_list(np_refi)=0
3809
                      end if
3810
                   end do
3811

    
3812
                   !!! Falta THE !!!
3813

    
3814
                case (6)
3815
                   !---- CEN ----!
3816
                   do j=1,3
3817
                      if (molec%lxcentre(j) ==0) then
3818
                         np_refi=np_refi+1
3819
                         V_Vec(np_refi)=molec%xcentre(j)
3820
                         V_name(np_refi)=trim(code_nam(12+j))//"Mol"
3821
                         molec%mxcentre(j)=1.0
3822
                         molec%lxcentre(j)=np_refi
3823
                         V_Bounds(1,np_refi)=xl
3824
                         V_Bounds(2,np_refi)=xu
3825
                         V_Bounds(3,np_refi)=xs
3826
                         V_BCon(np_refi)=ic
3827
                         V_list(np_refi)=0
3828
                      end if
3829
                   end do
3830

    
3831
                case (7)
3832
                   !---- ORI ----!
3833
                   do j=1,3
3834
                      if (molec%lOrient(j) ==0) then
3835
                         np_refi=np_refi+1
3836
                         V_Vec(np_refi)=molec%orient(j)
3837
                         V_name(np_refi)=trim(code_nam(15+j))//"Mol"
3838
                         molec%mOrient(j)=1.0
3839
                         molec%lOrient(j)=np_refi
3840
                         V_Bounds(1,np_refi)=xl
3841
                         V_Bounds(2,np_refi)=xu
3842
                         V_Bounds(3,np_refi)=xs
3843
                         V_BCon(np_refi)=ic
3844
                         V_list(np_refi)=0
3845
                      end if
3846
                   end do
3847

    
3848
                case (8)
3849
                   !---- THE ----!
3850

    
3851
                   !!! Not yet implemented !!!
3852
             end select
3853
       end select
3854

    
3855
       return
3856
    End Subroutine Fill_RefCodes_Molec
3857

    
3858
    !!--++
3859
    !!--++ Subroutine Fill_RefCodes_Magdom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Mag_dom)
3860
    !!--++    integer,                       intent(in)     :: Key
3861
    !!--++    character(len=*),              intent(in)     :: Dire
3862
    !!--++    integer,                       intent(in)     :: Na
3863
    !!--++    integer,                       intent(in)     :: Nb
3864
    !!--++    real(kind=cp),                 intent(in)     :: Xl
3865
    !!--++    real(kind=cp),                 intent(in)     :: Xu
3866
    !!--++    real(kind=cp),                 intent(in)     :: Xs
3867
    !!--++    integer,                       intent(in)     :: Ic
3868
    !!--++    type(Magnetic_Domain_type),    intent(in out) :: Mag_dom
3869
    !!--++
3870
    !!--++ Write on Vectors the Information for Magnetic Domains
3871
    !!--++ Created: February - 2012
3872
    !!--++
3873
    Subroutine Fill_RefCodes_Magdom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Mag_dom)
3874
       !---- Arguments ----!
3875
       integer,                       intent(in)     :: Key
3876
       character(len=*),              intent(in)     :: Dire
3877
       integer,                       intent(in)     :: Na
3878
       integer,                       intent(in)     :: Nb
3879
       real(kind=cp),                 intent(in)     :: Xl
3880
       real(kind=cp),                 intent(in)     :: Xu
3881
       real(kind=cp),                 intent(in)     :: Xs
3882
       integer,                       intent(in)     :: Ic
3883
       type(Magnetic_Domain_type),    intent(in out) :: Mag_dom
3884

    
3885
       !---- Local variables ----!
3886
       integer           :: i,j,nc
3887

    
3888
       call init_err_refcodes()
3889
       if (Na <= 0) then
3890
          err_refcodes=.true.
3891
          ERR_RefCodes_Mess="Number of domains no defined"
3892
          return
3893
       end if
3894

    
3895
       select case (dire)
3896
          !---- FIX Directive ----!
3897
          case ("fix")
3898

    
3899
             select case (key)
3900
                case (0)
3901
                   err_refcodes=.true.
3902
                   ERR_RefCodes_Mess="Option not defined"
3903
                   return
3904

    
3905
                case (4)
3906
                   !---- Magdomain ----!
3907
                      if (Mag_Dom%Lpop(nb,na) /=0) then
3908
                         nc=Mag_Dom%Lpop(nb,na)
3909
                         call Delete_RefCodes(nc,Mag_dom)
3910
                      end if
3911

    
3912
                case (2:3)
3913
                   err_refcodes=.true.
3914
                   ERR_RefCodes_Mess="Option not defined"
3915
                   return
3916

    
3917
                case (5:)
3918
                   err_refcodes=.true.
3919
                   ERR_RefCodes_Mess="Option not defined"
3920
                   return
3921
             end select
3922

    
3923
          !---- VARY Directive ----!
3924
          case ("var")
3925

    
3926
              select case (key)
3927
                case (0)
3928
                   err_refcodes=.true.
3929
                   ERR_RefCodes_Mess="Option not defined"
3930
                   return
3931

    
3932
                case (4)
3933
                   !---- Magd ----!
3934
                      if (Mag_Dom%Lpop(nb,na) ==0) then
3935
                            np_refi=np_refi+1
3936
                            Mag_Dom%Lpop(nb,na)=np_refi
3937
                            Mag_Dom%Mpop(nb,na)=1.0
3938

    
3939
                         if (Mag_Dom%Lpop(nb,na) == np_refi) then
3940
                            V_Vec(np_refi)=Mag_Dom%pop(nb,na)
3941
                            V_Name(np_refi)=trim(Mag_dom%lab(nb,na))
3942
                            V_Bounds(1,np_refi)=xl
3943
                            V_Bounds(2,np_refi)=xu
3944
                            V_Bounds(3,np_refi)=xs
3945
                            V_BCon(np_refi)=ic
3946
                            V_list(np_refi)=0
3947
                         else
3948
                            np_refi=np_refi-1
3949
                         end if
3950
                      end if
3951

    
3952
                case (2:3)
3953
                   err_refcodes=.true.
3954
                   ERR_RefCodes_Mess="Option not defined"
3955
                   return
3956

    
3957
                case (5:)
3958
                   err_refcodes=.true.
3959
                   ERR_RefCodes_Mess="Option not defined"
3960
                   return
3961
              end select
3962
       end select
3963

    
3964
       return
3965
    End Subroutine Fill_RefCodes_Magdom
3966

    
3967
    !!--++
3968
    !!--++  Subroutine Get_Atombet_Ctr(X,Betas,Spgr,Codini,Icodes,Multip,Ord,Ss,Ipr)
3969
    !!--++     real(kind=cp), dimension(3),             intent(in    ) :: X         !Atom position (fractional coordinates)
3970
    !!--++     real(kind=cp), dimension(6),             intent(in out) :: Betas     !Anisotropic temperature factors
3971
    !!--++     type(Space_Group_type),                  intent(in    ) :: Spgr      !Space Group
3972
    !!--++     Integer,                                 intent(in out) :: Codini    !Last attributed parameter
3973
    !!--++     Integer, dimension(6),                   intent(in out) :: Icodes    !codewords for betas only number
3974
    !!--++     real(kind=cp), dimension(6),             intent(in out) :: Multip    !Multipliers
3975
    !!--++     integer,                       optional, intent(in    ) :: Ord       !Order of the stabilizer
3976
    !!--++     integer, dimension(:),         optional, intent(in    ) :: Ss        !Pointer to SymmOp. of stabilizer
3977
    !!--++     integer,                       optional, intent(in    ) :: Ipr       !Printing unit for debug
3978
    !!--++
3979
    !!--++  Subroutine to get the appropriate constraints in the refinement codes of
3980
    !!--++  anisotropic atomic displacement(thermal) parameters.
3981
    !!--++  New algorithm based in the Wigner theorem.
3982
    !!--++  The matrix Bet = Sum { R Beta RT} displays the symmetry constraints to be
3983
    !!--++  applied to the anisotropic temperature factors. The sum runs over all rotational
3984
    !!--++  symmetry operators of the stabilizer of the particular atom position in the given
3985
    !!--++  space group.
3986
    !!--++  There are a total of 29 kind of relations that may appear in the Bet matrix:
3987
    !!--++
3988
    !!--++     1    A A A 0   0   0  -> m-3m, -43m, 432, m-3,23, 3[111].2[001]
3989
    !!--++     2    A A C 0   0   0  -> 4/mmm, -42m, 4mm, 422, 4/m, -4,4, 4[001]
3990
    !!--++     3    A B A 0   0   0  -> 4[010]
3991
    !!--++     4    A B B 0   0   0  -> 4[100]
3992
    !!--++     5    A A A D   D   D  -> -3m, 3m, 32, -3, 3   3[111]
3993
    !!--++     6    A A A D  -D  -D  -> 3[11-1]
3994
    !!--++     7    A A A D  -D   D  -> 3[1-11]
3995
    !!--++     8    A A A D   D  -D  -> 3[-111]
3996
    !!--++     9    A A C A/2 0   0  -> 6/mmm, -6m2, 6mm, 622, 6/m, 6,-6,-3m, 32,-3, 3:  h 3[001]
3997
    !!--++    10    A B C 0   0   0  -> mmm, mm2, 222  2[001] 2[100]
3998
    !!--++    11    A A C D   0   0  -> 2[001], 2[110]    w
3999
    !!--++    12    A B A 0   E   0  -> 2[010], 2[101]
4000
    !!--++    13    A B B 0   0   F  -> 2[100], 2[011]
4001
    !!--++    14    A B C B/2 0   0  -> 2[001], 2[100]    h
4002
    !!--++    15    A B C A/2 0   0  -> 2[001], 2[010]    h
4003
    !!--++    16    A B C D   0   0  -> 2/m, m, 2: 2[001] w
4004
    !!--++    17    A B C 0   E   0  -> 2[010]
4005
    !!--++    18    A B C 0   0   F  -> 2[100]
4006
    !!--++    19    A A C D   E  -E  -> 2[110]            w
4007
    !!--++    20    A A C D   E   E  -> 2[1-10]           w
4008
    !!--++    21    A B A D   E  -D  -> 2[101]
4009
    !!--++    22    A B A D   E   D  -> 2[10-1]
4010
    !!--++    23    A B B D  -D   F  -> 2[011]
4011
    !!--++    24    A B B D   D   F  -> 2[01-1]
4012
    !!--++    25    A B C B/2 F/2 F  -> 2[100]            h
4013
    !!--++    26    A B C A/2 0   F  -> 2[210]            h
4014
    !!--++    27    A B C B/2 E   0  -> 2[120]            h
4015
    !!--++    28    A B C A/2 E   E/2-> 2[010]            h
4016
    !!--++    29    A B C D   E   F  -> 1, -1
4017
    !!--++
4018
    !!--++   Updated: 14 July 2011
4019
    !!--++
4020
    !!
4021

    
4022
    Subroutine Get_Atombet_Ctr(X,Betas,Spgr,Codini,Icodes,Multip,Ord,Ss,Ipr)
4023
       real(kind=cp), dimension(3),             intent(in    ) :: X
4024
       real(kind=cp), dimension(6),             intent(in out) :: Betas
4025
       type(Space_Group_type),                  intent(in    ) :: Spgr
4026
       Integer,                                 intent(in out) :: Codini
4027
       Integer, dimension(6),                   intent(in out) :: Icodes
4028
       real(kind=cp), dimension(6),             intent(in out) :: Multip
4029
       integer,                       optional, intent(in    ) :: Ord
4030
       integer, dimension(:),         optional, intent(in    ) :: Ss
4031
       integer,                       optional, intent(in    ) :: Ipr
4032

    
4033

    
4034
       ! Local variables
4035
       character (len=1), dimension(6)   :: cdd
4036
       integer                           :: i,j,order
4037
       integer,           dimension(48)  :: ss_ptr
4038
       integer,           dimension(6)   :: codd
4039
       integer,           dimension(3,3) :: Rsym
4040
       real(kind=cp),     dimension(3,3) :: bet,bett,Rs
4041
       real(kind=cp),     dimension(6)   :: cod
4042
       real(kind=cp),     dimension(3,48):: atr
4043
       real(kind=cp),     parameter      :: epss=0.01_cp
4044

    
4045
       cod=real(icodes)
4046

    
4047
       do j=1,6
4048
          if (cod(j) < 1.0 .and. abs(multip(j)) > epss)  then
4049
             codini=codini+1
4050
             cod(j) = real(codini)
4051
          end if
4052
       end do
4053

    
4054
       if (present(ord) .and. present(ss)) then
4055
          order=ord
4056
          ss_ptr(1:order) = ss(1:ord)
4057
       else
4058
          call get_stabilizer(x,Spgr,order,ss_ptr,atr)
4059
       end if
4060

    
4061
       bet=reshape((/17.0, 7.0,3.0,  &
4062
                     7.0,13.0,5.0,  &
4063
                     3.0, 5.0,11.0/),(/3,3/))
4064
       bett=bet
4065
       if (order > 1 ) then
4066
          do j=2,order
4067
             Rsym=Spgr%SymOp(ss_ptr(j))%Rot
4068
             Rs=real(Rsym)
4069
             bett=bett+ matmul(Rs,matmul(bet,transpose(Rs)))
4070
          end do
4071
       end if
4072
       Rsym=nint(1000.0*bett)
4073
       codd=(/Rsym(1,1),Rsym(2,2),Rsym(3,3),Rsym(1,2),Rsym(1,3),Rsym(2,3)/)
4074
       cdd=(/'a','b','c','d','e','f'/)
4075
       multip=1.0
4076
       !Search systematically all the possible constraints
4077

    
4078
       if(codd(1) == codd(2) .and. codd(1) == codd(3)) then ! a a a
4079
         if(codd(4) == codd(5) .and. codd(4) == codd(6) ) then ! a a a d d d
4080
           if(codd(4) == 0) then
4081
             cdd=(/'a','a','a','0','0','0'/)     ! 1 A A A 0   0   0
4082
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4083
             betas(4:6)=0.0
4084
             betas(2:3)=betas(1)
4085
             cod(2:3)=cod(1); cod(4:6)=0.0
4086
           else
4087
             cdd=(/'a','a','a','d','d','d'/)     ! 5 A A A D   D   D
4088
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4089
             betas(5:6)=betas(4)
4090
             betas(2:3)=betas(1)
4091
             cod(2:3)=cod(1); cod(5:6)=cod(4)
4092
           end if
4093
         else if(codd(4) == -codd(5) .and. codd(4) == -codd(6) ) then !a a a d -d -d
4094
           cdd=(/'a','a','a','d','d','d'/)       ! 6 A A A D  -D  -D
4095
           multip=(/1.0,1.0,1.0,1.0,-1.0,-1.0/)
4096
           betas(5:6)=-betas(4)
4097
           betas(2:3)=betas(1)
4098
           cod(2:3)=cod(1); cod(5:6)=cod(4)
4099
         else if(codd(4) == -codd(5) .and. codd(4) ==  codd(6) ) then !a a a d -d  d
4100
           cdd=(/'a','a','a','d','d','d'/)       ! 7 A A A D  -D   D
4101
           multip=(/1.0,1.0,1.0,1.0,-1.0, 1.0/)
4102
           betas(5)=-betas(4); betas(6)=betas(4)
4103
           betas(2:3)=betas(1)
4104
           cod(2:3)=cod(1); cod(5:6)= cod(4)
4105
         else if(codd(4) ==  codd(5) .and. codd(4) == -codd(6) ) then !a a a d  d -d
4106
           cdd=(/'a','a','a','d','d','d'/)       ! 8 A A A D   D  -D
4107
           multip=(/1.0,1.0,1.0,1.0, 1.0,-1.0/)
4108
           betas(6)=-betas(4); betas(5)=betas(4)
4109
           betas(2:3)=betas(1)
4110
           cod(2:3)=cod(1); cod(5:6)= cod(4)
4111
         end if
4112

    
4113
       else if(codd(1) == codd(2)) then ! a a c
4114
         if(codd(4) == codd(5) .and. codd(4) == codd(6) .and. codd(4) == 0) then ! a a c 0 0 0
4115
             cdd=(/'a','a','c','0','0','0'/)     ! 2 A A C 0   0   0
4116
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4117
             betas(4:6)=0.0
4118
             betas(2)=betas(1)
4119
             cod(2)=cod(1); cod(4:6)= 0.0
4120
         else if(codd(5) == codd(6) .and. codd(5) == 0) then ! a a c x 0 0
4121
             if(codd(4) == codd(1)/2) then
4122
               cdd=(/'a','a','c','a','0','0'/)     ! 9 A A C A/2 0   0
4123
               multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
4124
               betas(5:6)=0.0; betas(4)=betas(1)*0.5
4125
               betas(2)=betas(1)
4126
               cod(2)=cod(1); cod(4)= cod(1); cod(5:6)=0.0
4127
             else
4128
               cdd=(/'a','a','c','d','0','0'/)     !11 A A C D   0   0
4129
               multip=(/1.0,1.0,1.0,1.0,0.0,0.0/)
4130
               betas(5:6)=0.0
4131
               betas(2)=betas(1)
4132
               cod(2)=cod(1); cod(5:6)=0.0
4133
             end if
4134
         else
4135
             if(codd(5) == codd(6)) then  ! a a c d e e
4136
               cdd=(/'a','a','c','d','e','e'/)     !20 A A C D   E   E
4137
               multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4138
               betas(6)=betas(5)
4139
               betas(2)=betas(1)
4140
               cod(2)=cod(1); cod(6)=cod(5)
4141
             else if(codd(5) == -codd(6)) then  ! a a c d e -e
4142
               cdd=(/'a','a','c','d','e','e'/)     !19 A A C D   E  -E
4143
               multip=(/1.0,1.0,1.0,1.0,1.0,-1.0/)
4144
               betas(6)=-betas(5)
4145
               betas(2)=betas(1)
4146
               cod(2)=cod(1); cod(6)=cod(5)
4147
             end if
4148
         end if
4149

    
4150
       else if(codd(1) == codd(3)) then ! a b a
4151
         if(codd(4) == codd(6)) then    ! a b a d x d
4152
           if(codd(4) == 0) then  ! a b a 0 x 0
4153
             if(codd(5) == 0) then ! a b a 0 0 0
4154
               cdd=(/'a','b','a','0','0','0'/)     ! 3 A B A 0   0   0
4155
               multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4156
               betas(4:6)=0.0
4157
               betas(3)=betas(1)
4158
               cod(3)=cod(1); cod(4:6)=0.0
4159
             else                  ! a b a 0 e 0
4160
               cdd=(/'a','b','a','0','e','0'/)     !12 A B A 0   E   0
4161
               multip=(/1.0,1.0,1.0,0.0,1.0,0.0/)
4162
               betas(4)=0.0;  betas(6)=0.0
4163
               betas(3)=betas(1)
4164
               cod(3)=cod(1); cod(4)=0.0;  cod(6)=0.0
4165
             endif
4166
           else  !! a b a d e d
4167
             cdd=(/'a','b','a','d','e','d'/)       !22 A B A D   E   D
4168
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4169
             betas(6)=betas(4)
4170
             betas(3)=betas(1)
4171
             cod(3)=cod(1); cod(6)=cod(4)
4172
          end if
4173

    
4174
         else if(codd(4) == -codd(6)) then ! a b a d e -d
4175
           cdd=(/'a','b','a','d','e','d'/)         !21 A B A D   E  -D
4176
           multip=(/1.0,1.0,1.0,1.0,1.0,-1.0/)
4177
           betas(6)=-betas(4)
4178
           betas(3)=betas(1)
4179
           cod(3)=cod(1); cod(6)=cod(4)
4180
         end if
4181

    
4182
       else if(codd(2) == codd(3)) then ! a b b
4183
         if(codd(4) == codd(5)) then    ! a b b d d x
4184
           if(codd(4) == 0) then  ! a b b 0 0 x
4185
             if(codd(6) == 0) then ! a b b 0 0 0
4186
               cdd=(/'a','b','b','0','0','0'/)     ! 4 A B B 0   0   0
4187
               multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4188
               betas(4:6)=0.0
4189
               betas(3)=betas(2)
4190
               cod(3)=cod(2); cod(4:6)=0.0
4191
             else                  ! a b b 0 0 f
4192
               cdd=(/'a','b','b','0','0','f'/)     !13 A B B 0   0   F
4193
               multip=(/1.0,1.0,1.0,0.0,0.0,1.0/)
4194
               betas(4:5)=0.0
4195
               betas(3)=betas(2)
4196
               cod(3)=cod(2); cod(4:5)=0.0
4197
             endif
4198
           else  !! a b b d d f
4199
             cdd=(/'a','b','b','d','d','f'/)       !24 A B B D   D   F
4200
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4201
             betas(5)=betas(4)
4202
             betas(3)=betas(2)
4203
             cod(3)=cod(2); cod(5)=cod(4)
4204
           end if
4205
         else if(codd(4) == -codd(5)) then ! a b b d -d e
4206
           cdd=(/'a','b','b','d','d','f'/)         !23 A B B D  -D   F
4207
           multip=(/1.0,1.0,1.0,1.0,-1.0,1.0/)
4208
           betas(5)=-betas(4)
4209
           betas(3)=betas(2)
4210
           cod(3)=cod(2); cod(5)=cod(4)
4211
         end if
4212

    
4213
       else !Now a /= b /= c
4214

    
4215
         if(codd(4) == codd(5) .and. codd(4) == 0) then ! a b c 0 0 x
4216
           if(codd(6) == 0) then ! a b c 0 0 0
4217
             cdd=(/'a','b','c','0','0','0'/)          !10 A B C 0   0   0
4218
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4219
             betas(4:6)=0.0
4220
             cod(4:6)=0.0
4221
           else
4222
             cdd=(/'a','b','c','0','0','f'/)          !18 A B C 0   0   F
4223
             multip=(/1.0,1.0,1.0,0.0,0.0,1.0/)
4224
             betas(4:5)=0.0
4225
             cod(4:5)=0.0
4226
           end  if
4227
         else if(codd(5) == codd(6) .and. codd(5) == 0) then  ! a b c x 0 0
4228
           if(codd(4) == codd(1)/2) then ! a b c a/2 0 0
4229
             cdd=(/'a','b','c','a','0','0'/)          !15 A B C A/2 0   0
4230
             multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
4231
             betas(5:6)=0.0; betas(4)=betas(1)*0.5
4232
             cod(4)=cod(1); cod(5:6)=0.0
4233
           else if(codd(4) == codd(2)/2) then    !a b c b/2 0 0
4234
             cdd=(/'a','b','c','b','0','0'/)          !14 A B C B/2 0   0
4235
             multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
4236
             betas(5:6)=0.0; betas(4)=betas(2)*0.5
4237
             cod(4)=cod(2); cod(5:6)=0.0
4238
           else
4239
             cdd=(/'a','b','c','d','0','0'/)          !16 A B C D   0   0
4240
             multip=(/1.0,1.0,1.0,1.0,0.0,0.0/)
4241
             betas(5:6)=0.0
4242
             cod(5:6)=0.0
4243
           end  if
4244
         else if(codd(4) == codd(6) .and. codd(4) == 0) then !a b c 0 e 0
4245
           cdd=(/'a','b','c','0','e','0'/)            !17 A B C 0   E   0
4246
           multip=(/1.0,1.0,1.0,0.0,1.0,0.0/)
4247
           betas(4)=0.0; betas(6)=0.0
4248
           cod(4)=0.0; cod(6)=0.0
4249
         else if(codd(4) == codd(1)/2 .and. codd(5) == 0) then !a b c a/2 0 f
4250
           cdd=(/'a','b','c','a','0','f'/)            !26 A B C A/2 0   F
4251
           multip=(/1.0,1.0,1.0,0.5,0.0,1.0/)
4252
           betas(4)=betas(1)*0.5; betas(5)=0.0
4253
           cod(4)=cod(1); cod(5)=0.0
4254
         else if(codd(4) == codd(2)/2 .and. codd(6) == 0) then !a b c b/2 e 0
4255
           cdd=(/'a','b','c','b','e','0'/)            !27 A B C B/2 E   0
4256
           multip=(/1.0,1.0,1.0,0.5,1.0,0.0/)
4257
           betas(4)=betas(2)*0.5; betas(6)=0.0
4258
           cod(4)=cod(2); cod(6)=0.0
4259
         else if(codd(4) == codd(2)/2 .and. codd(5) == codd(6)/2) then !a b c b/2 f/2 f
4260
           cdd=(/'a','b','c','b','f','f'/)            !25 A B C B/2 F/2 F
4261
           multip=(/1.0,1.0,1.0,0.5,0.5,1.0/)
4262
           betas(4)=betas(2)*0.5; betas(5)=betas(6)*0.5
4263
           cod(4)=cod(2); cod(5)=cod(6)
4264
         else if(codd(4) == codd(1)/2 .and. codd(6) == codd(5)/2) then !a b c a/2 e e/2
4265
           cdd=(/'a','b','c','a','e','e'/)            !28 A B C A/2 E   E/2
4266
           multip=(/1.0,1.0,1.0,0.5,1.0,0.5/)
4267
           betas(4)=betas(1)*0.5; betas(6)=betas(5)*0.5
4268
           cod(4)=cod(1); cod(6)=cod(5)
4269
         else
4270
           cdd=(/'a','b','c','d','e','f'/)            !29 A B C D   E   F
4271
           multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4272
         end if
4273
       end if
4274

    
4275
       do j=1,6
4276
          if (multip(j) < epss .or. cdd(j) == "0" ) then
4277
             icodes(j) = 0
4278
          else
4279
             icodes(j) = nint(cod(j))
4280
          end if
4281
       end do
4282

    
4283
       if(present(Ipr)) then
4284
         Write(Ipr,'(a,6i5)')           '     Codes on Betas       :  ',Icodes
4285
         Write(Ipr,'(a,6(a,1x),6f7.3)') '     Codes and multipliers:  ',cdd,multip
4286
         Write(Ipr,'(a)')               '     Beta_TOT matrix:  '
4287
         Do I=1,3
4288
          Write(Ipr,'(a,3f12.4)')       '                      ',bett(i,:)
4289
         End Do
4290
       end if
4291
       return
4292
    End Subroutine Get_Atombet_Ctr
4293

    
4294
    !!--++
4295
    !!--++  Subroutine Get_Atompos_Ctr(X,Spgr,Codini,Icodes,Multip,Ord,Ss,Att,Ipr)
4296
    !!--++     real(kind=cp), dimension(3),       intent(in    ) :: x      !Atom position (fractional coordinates)
4297
    !!--++     type(Space_Group_type),            intent(in    ) :: Spgr   !Space Group
4298
    !!--++     Integer,                           intent(in out) :: codini !Last attributed parameter
4299
    !!--++     integer,       dimension(3),       intent(in out) :: Icodes
4300
    !!--++     real(kind=cp), dimension(3),       intent(in out) :: Multip
4301
    !!--++     integer,                 optional, intent(in    ) :: Ord
4302
    !!--++     integer, dimension(:),   optional, intent(in    ) :: Ss
4303
    !!--++     integer, dimension(:,:), optional, intent(in    ) :: Atr
4304
    !!--++     integer,                 optional, intent(in    ) :: Ipr
4305
    !!--++
4306
    !!--++     (Private)
4307
    !!--++     Subroutine to get the appropriate constraints in the refinement codes of
4308
    !!--++     atoms positions. The algorithm is based in an analysis of the symbol generated
4309
    !!--++     for the symmetry elements of the operators belonging to the stabilizer of the
4310
    !!--++     atom position. This routine operates with splitted codes in the sense of
4311
    !!--++     FullProf rules
4312
    !!--++
4313
    !!--++     Updated: July - 2011 (JRC, the old subroutine has been completely changed)
4314
    !!
4315
    Subroutine Get_Atompos_Ctr(X,Spgr,Codini,ICodes,Multip,Ord,Ss,Att,Ipr)
4316
       real(kind=cp), dimension(3),            intent(in)     :: X
4317
       type(Space_Group_type),                 intent(in)     :: Spgr
4318
       Integer,                                intent(in out) :: Codini
4319
       Integer,       dimension(3),            intent(in out) :: ICodes
4320
       real(kind=cp), dimension(3),            intent(in out) :: Multip
4321
       integer,                       optional,intent(in)     :: Ord
4322
       integer, dimension(:),         optional,intent(in)     :: Ss
4323
       real(kind=cp), dimension(:,:), optional,intent(in)     :: Att
4324
       integer,                       optional,intent(in)     :: Ipr
4325

    
4326
       ! Local variables
4327
       integer                          :: i,j,k,order,L,L1,L2,ipar,j1
4328
       integer,          dimension(3,3) :: RSym
4329
       integer,          dimension(48)  :: ss_ptr
4330
       real(kind=cp),    dimension(3,48):: atr
4331
       real(kind=cp),    dimension(3)   :: tr
4332

    
4333
       character(len=40)                :: symbol,tsymbol,sym_symb
4334
       Character(len=10), dimension(3)  :: nsymb
4335
       Character(len=3),  dimension(3)  :: ssymb
4336
       real(kind=cp),     parameter     :: epss=0.001
4337

    
4338
       if(present(ord) .and. present(ss) .and. present(att)) then
4339
         order=ord
4340
         ss_ptr(1:order) = ss(1:ord)
4341
         atr(:,1:order)  = att(:,1:ord)
4342
       else
4343
         call get_stabilizer(x,Spgr,order,ss_ptr,atr)
4344
       end if
4345

    
4346
       !If codes were not assigned with explicit number
4347
       !attribute numbers bigger than initial Codini
4348

    
4349
       do j=1,3
4350
          if(Icodes(j) < 1  .and. abs(multip(j)) > epss)  then
4351
               codini = codini+1
4352
               Icodes(j) = codini
4353
          end if
4354
       end do
4355

    
4356
       ssymb=(/"  x","  y","  z"/)
4357

    
4358
       if(present(Ipr)) write(unit=Ipr,fmt='(/a,3f10.5)')  ' => Atom Position:',x
4359

    
4360
       if(order > 1 ) then  !A constraint in position must exist
4361

    
4362
          if(present(Ipr)) write(unit=Ipr,fmt='(a)')   ' => List of symmetry element of the stabilizer without identity:'
4363

    
4364
          do k=2,order
4365
           symbol=" "
4366
           Rsym=Spgr%SymOp(ss_ptr(k))%Rot
4367
           tr=Spgr%SymOp(ss_ptr(k))%tr + atr(:,k)
4368
           call Get_SymSymb(Rsym,tr,Sym_Symb)
4369
           call symmetry_symbol(Sym_Symb,tsymbol)
4370
             i=index(tsymbol,";")
4371
             if(i /= 0) then
4372
               symbol=tsymbol(1:i-1)
4373
               call Read_Xsym(tsymbol(i+1:),1,Rsym,Tr,.false.)
4374
               if(sum(abs(x-tr)) < epss) then
4375
                  ssymb=(/"  0","  0","  0"/)
4376
                  if(present(Ipr)) then
4377
                    write(unit=Ipr,fmt="(a,i2,a,t20,a,t55,a,t90,4a)") "     Operator ",k,": ", &
4378
                    trim(Sym_Symb),trim(tsymbol),"  ssymb:" ,(ssymb(j)//"  ",j=1,3)
4379
                  end if
4380
                  cycle
4381
               end if
4382
             else
4383
               symbol=tsymbol
4384
             end if
4385
             ipar=index(symbol,")")              !Translation element appears before position
4386
             L =index(symbol(ipar+1:)," ")+ipar  !Position of the first blank after translation
4387
             L1=index(symbol(ipar+1:),",")+ipar  !Position of the first comma after translation
4388
             L2=index(symbol(L1+1:),",")+L1      !Position of the second comma
4389
             if(L1 == 0) L1=1
4390
             if(L2 == 0) L2=1
4391
             if(L  == 0) L=1
4392

    
4393
             !Construct a new symbol that estabish automatically the constraints
4394
             nsymb = (/symbol(L+1:L1-1),symbol(L1+1:L2-1),symbol(L2+1:)/)
4395

    
4396
             do i=1,3
4397
                 do j=1,10  !Delete unwanted symbols (keep only x,y,z,2 and -
4398
                    if(nsymb(i)(j:j) == " ") cycle
4399
                    if(nsymb(i)(j:j) /= "x" .and. nsymb(i)(j:j) /= "y" .and. &
4400
                       nsymb(i)(j:j) /= "z" .and. nsymb(i)(j:j) /= "-" .and. &
4401
                       nsymb(i)(j:j) /= "2" ) nsymb(i)(j:j)=" "
4402
                 end do
4403
                 if(len_trim(nsymb(i))  == 0 .or. (index(nsymb(i),"x") == 0 .and. &
4404
                    index(nsymb(i),"y") == 0 .and. index(nsymb(i),"z") == 0  ) ) then
4405
                    ssymb(i)="  0"
4406
                    cycle
4407
                 end if
4408
                 !Now remove 2s on the right of x,y, or z
4409
                 j1=index(nsymb(i),"2")
4410
                 if( j1 /= 0) then
4411
                    if(len_trim(nsymb(i)) == j1) nsymb(i)=nsymb(i)(1:j1-1)
4412
                 end if
4413
                 !Now remove -s on the right of x,y, or z
4414
                 j1=index(nsymb(i),"-")
4415
                 if( j1 /= 0) then
4416
                    if(len_trim(nsymb(i)) == j1) nsymb(i)=nsymb(i)(1:j1-1)
4417
                 end if
4418
                 nsymb(i)= adjustl(nsymb(i))
4419
             end do
4420

    
4421
             if(ssymb(1) /= "  0" .and. ssymb(1) /= "  a") then
4422
                ssymb(1)= nsymb(1)
4423
                ssymb(1)= adjustr(ssymb(1))
4424
             end if
4425

    
4426
             if(ssymb(2) /= "  0" .and. ssymb(2) /= "  a" .and. ssymb(2) /= "  b" .and. &
4427
                ssymb(2) /= " -a" .and. ssymb(2) /= " 2a"   ) then
4428
                ssymb(2) = nsymb(2)
4429
                ssymb(2) = adjustr(ssymb(2))
4430
             end if
4431

    
4432
             if(ssymb(3) /= "  0" .and. ssymb(3) /= "  a" .and. ssymb(3) /= "  b" .and. &
4433
                ssymb(3) /= "  c" .and. ssymb(3) /= " 2a" .and. ssymb(3) /= " 2b" .and. &
4434
                ssymb(3) /= " -a" .and. ssymb(3) /= " -b") then
4435
                ssymb(3) = nsymb(3)
4436
                ssymb(3) = adjustr(ssymb(3))
4437
             end if
4438

    
4439
             do i=1,3
4440
                if(ssymb(i)(3:3) == "x")  ssymb(i)(3:3) = "a"
4441
             end do
4442
             do i=1,3
4443
                if(ssymb(i)(3:3) == "y")  ssymb(i)(3:3) = "b"
4444
             end do
4445
             do i=1,3
4446
                if(ssymb(i)(3:3) == "z")  ssymb(i)(3:3) = "c"
4447
             end do
4448
             if(present(Ipr)) then
4449
                write(unit=Ipr,fmt="(a,i2,a,t20,a,t55,a,t90,4a)") "     Operator ",k,": ", &
4450
                trim(Sym_Symb),trim(tsymbol),"  Ssymb:" ,(ssymb(j)//"  ",j=1,3)
4451
             end if
4452

    
4453
          end do !do k=1,order  over operators of the stabilizer
4454

    
4455
       else
4456
         ssymb=(/"  a","  b","  c"/)
4457

    
4458
       end if  !order > 1
4459

    
4460
       do i=1,3                  !Fixing codes
4461
         if(ssymb(i)=="  0") then
4462
           Icodes(i)=0
4463
           multip(i)=0.0
4464
         end if
4465
       end do
4466

    
4467
       if(index(ssymb(1),"a") /= 0) then
4468

    
4469
         do i=2,3  !Fixing codes
4470
           if(index(ssymb(i),"-a") /= 0) then
4471
             Icodes(i)=Icodes(1)
4472
             multip(i)=-multip(1)
4473
           else if(index(ssymb(i),"a") /= 0) then
4474
             Icodes(i)=Icodes(1)
4475
             multip(i)=multip(1)
4476

    
4477
             if(index(ssymb(i),"2") /= 0) then
4478
               multip(i)=2.0* multip(1)
4479
             else if(index(ssymb(1),"2") /= 0) then
4480
               multip(i)=0.5* multip(1)
4481
             end if
4482

    
4483
           end if
4484
         end do
4485
       else  !the x-coordinate is fixed, analyse y and z
4486
         if(index(ssymb(2),"b") /= 0 .and. index(ssymb(3),"b") /= 0) then
4487
           Icodes(3)=Icodes(2)
4488
           if(ssymb(2) == ssymb(3)) then
4489
             multip(3)= multip(2)
4490
           else if(ssymb(3) == " -b" .and. ssymb(2) == "  b") then
4491
             multip(3)= -multip(2)
4492
           else if(ssymb(3) == "  b" .and. ssymb(2) == " -b") then
4493
             multip(3)= -multip(2)
4494
           end if
4495
         end if
4496

    
4497
       end if !if(index(ssymb(1),"a") /= 0)
4498

    
4499
       do j=1,3
4500
         if(abs(multip(j)) < epss) then
4501
           Icodes(j) = 0
4502
         end if
4503
       end do
4504
       if(present(Ipr)) then
4505
         write(unit=Ipr,fmt="(a,3i5)")    "     Codes positions: ",Icodes
4506
         write(unit=Ipr,fmt="(a,3f5.1)")  "     Multipliers    : ",multip
4507
         write(unit=Ipr,fmt="(5a)")       "     Codes   string : ( ",(ssymb(j),j=1,3) ," )"
4508
       end if
4509
       return
4510
    End Subroutine Get_Atompos_Ctr
4511

    
4512
    !!--++
4513
    !!--++ Subroutine Get_ConCodes_Line(Line,FAtom/FmAtom/MolCrys/Molec/MagDom)
4514
    !!--++    character(len=*),             intent(in)     :: Line
4515
    !!--++    integer,                      intent(in)     :: Nat
4516
    !!--++    type(Atom_List_Type),         intent(in out) :: FAtom
4517
    !!--++    or
4518
    !!--++    type(mAtom_List_Type),        intent(in out) :: FmAtom
4519
    !!--++    or
4520
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
4521
    !!--++    or
4522
    !!--++    type(molecule_type),          intent(in out) :: Molec
4523
    !!--++    or
4524
    !!--++    type(Magnetic_Domain_type),   intent(in out) :: Mag_Dom
4525
    !!--++
4526
    !!--++    (Private)
4527
    !!--++    Get the Constraints relations
4528
    !!--++
4529
    !!--++ Update: March - 2005
4530
    !!
4531

    
4532
    !!--++
4533
    !!--++ Subroutine Get_ConCodes_Line_FAtom(Line,FAtom)
4534
    !!--++    character(len=*),         intent(in)     :: Line
4535
    !!--++    integer,                  intent(in)     :: Nat
4536
    !!--++    type(Atom_List_Type),     intent(in out) :: FAtom
4537
    !!--++
4538
    !!--++ Overloaded
4539
    !!--++ Get the Constraints relations
4540
    !!--++
4541
    !!--++ Update: March - 2005
4542
    !!
4543
    Subroutine Get_ConCodes_Line_FAtom(Line,FAtom)
4544
       !---- Arguments ----!
4545
       character(len=*),     intent(in)     :: Line
4546
       type(Atom_List_Type), intent(in out) :: FAtom
4547

    
4548
       !---- Local variables ----!
4549
       character(len=20), dimension(30) :: label
4550
       integer                          :: j,ic,n,na,nb,nc,nd,npos
4551
       integer                          :: nl,nl2,iv
4552
       integer, dimension(1)            :: ivet
4553
       real(kind=cp)                    :: fac_0,fac_1
4554
       real(kind=cp),dimension(1)       :: vet
4555

    
4556
       call init_err_refcodes()
4557

    
4558
       nl=0
4559
       call getword(line,label,ic)
4560
       if (ic < 2) then
4561
          err_refcodes=.true.
4562
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
4563
          return
4564
       end if
4565

    
4566
       !---- Set Futher Information----!
4567
       !---- Na is the number of atom on List
4568
       !---- Nb is the key (X,Y,Z,Occ,...)
4569
       !---- Fac0 is the multiplier
4570
       !---- Nl is the number of refinement parameter
4571

    
4572
       npos=index(label(1),"_")
4573
       if (npos ==0) then
4574
          err_refcodes=.true.
4575
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4576
          return
4577
       end if
4578

    
4579
       na=0
4580
       do j=1,FAtom%Natoms
4581
          if (u_case(FAtom%atom(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4582
             na=j
4583
             exit
4584
          end if
4585
       end do
4586
       if (na == 0) then
4587
          err_refcodes=.true.
4588
          ERR_RefCodes_Mess=" Atom label not found for "//trim(line)
4589
          return
4590
       end if
4591

    
4592
       nb=0
4593
       do j=1,ncode
4594
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
4595
             nb=j
4596
             exit
4597
          end if
4598
       end do
4599
       if (nb == 0) then
4600
          err_refcodes=.true.
4601
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
4602
          return
4603
       end if
4604

    
4605
       select case (nb)
4606
          case ( 1:3)
4607
             fac_0=FAtom%atom(na)%mx(nb)
4608
                nl=FAtom%atom(na)%lx(nb)
4609
          case ( 4)
4610
             fac_0=FAtom%atom(na)%mbiso
4611
                nl=FAtom%atom(na)%lbiso
4612
          case ( 5)
4613
             fac_0=FAtom%atom(na)%mocc
4614
                nl=FAtom%atom(na)%locc
4615
          case ( 6:11)
4616
             fac_0=FAtom%atom(na)%mu(nb-5)
4617
                nl=FAtom%atom(na)%lu(nb-5)
4618
          case (12)
4619
             fac_0=FAtom%atom(na)%mu(1)
4620
          case (13:)
4621
             err_refcodes=.true.
4622
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4623
             return
4624
       end select ! nb
4625

    
4626
       if (nb < ncode) then
4627
          if (nl == 0) then
4628
             err_refcodes=.true.
4629
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
4630
             return
4631
          end if
4632
       end if
4633

    
4634
       !---- Set the rest elements in Contsraints ----!
4635
       n=1
4636
       do
4637
          n=n+1
4638
          if (n > ic) exit
4639

    
4640
          npos=index(label(n),"_")
4641
          if (npos ==0) then
4642
             err_refcodes=.true.
4643
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
4644
             return
4645
          end if
4646

    
4647
          nc=0
4648
          do j=1,FAtom%Natoms
4649
             if (u_case(FAtom%atom(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
4650
                nc=j
4651
                exit
4652
             end if
4653
          end do
4654
          if (nc == 0) then
4655
             err_refcodes=.true.
4656
             ERR_RefCodes_Mess=" Atom label not found for "//trim(label(n))
4657
             return
4658
          end if
4659

    
4660
          nd=0
4661
          do j=1,ncode
4662
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
4663
                nd=j
4664
                exit
4665
             end if
4666
          end do
4667
          if (nd == 0) then
4668
             err_refcodes=.true.
4669
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
4670
             return
4671
          end if
4672

    
4673
          !---- Is there a new multiplier?
4674
          n=n+1
4675
          call getnum(label(n),vet,ivet,iv)
4676
          if (iv == 1) then
4677
             fac_1=vet(1)
4678
          else
4679
             fac_1=fac_0
4680
             n=n-1
4681
          end if
4682

    
4683
          select case (nd)
4684
             case ( 1:3)
4685
                nl2=FAtom%atom(nc)%lx(nd)
4686
                call Delete_refCodes(nl2,FAtom)
4687
                FAtom%atom(nc)%mx(nd)=fac_1
4688
                FAtom%atom(nc)%lx(nd)=nl
4689
             case ( 4)
4690
                nl2=FAtom%atom(nc)%lbiso
4691
                call Delete_refCodes(nl2,FAtom)
4692
                FAtom%atom(nc)%mbiso=fac_1
4693
                FAtom%atom(nc)%lbiso=nl
4694
             case ( 5)
4695
                nl2=FAtom%atom(nc)%locc
4696
                call Delete_refCodes(nl2,FAtom)
4697
                FAtom%atom(nc)%mocc=fac_1
4698
                FAtom%atom(nc)%locc=nl
4699
             case ( 6:11)
4700
                nl2=FAtom%atom(nc)%lu(nd-5)
4701
                call Delete_refCodes(nl2,FAtom)
4702
                FAtom%atom(nc)%mu(nd-5)=fac_1
4703
                FAtom%atom(nc)%lu(nd-5)=nl
4704
             case (12)
4705
                do j=1,6
4706
                   nl2=FAtom%atom(nc)%lu(j)
4707
                   call Delete_refCodes(nl2,FAtom)
4708
                   FAtom%atom(nc)%mu(j)=fac_1
4709
                   FAtom%atom(nc)%lu(j)=FAtom%atom(na)%lu(j)
4710
                end do
4711
                np_cons=np_cons+5
4712
             case (13:)
4713
                err_refcodes=.true.
4714
                ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4715
                return
4716
          end select ! nb
4717

    
4718
          np_cons=np_cons+1
4719

    
4720
       end do
4721

    
4722
       return
4723
    End Subroutine Get_ConCodes_Line_FAtom
4724

    
4725
    !!--++
4726
    !!--++ Subroutine Get_ConCodes_Line_FmAtom(Line,FmAtom)
4727
    !!--++    character(len=*),         intent(in)     :: Line
4728
    !!--++    integer,                  intent(in)     :: Nat
4729
    !!--++    type(mAtom_List_Type),    intent(in out) :: FmAtom
4730
    !!--++
4731
    !!--++ Get the Magnetic Constraints Relations: Presently only 'equal'
4732
    !!--++ mag clone of Get_ConCodes_Line_FAtom
4733
    !!--++ Created: December - 2011
4734
    !!
4735
    Subroutine Get_ConCodes_Line_FmAtom(Line,FmAtom)
4736
       !---- Arguments ----!
4737
       character(len=*),     intent(in)     :: Line
4738
       type(mAtom_List_Type),intent(in out) :: FmAtom
4739

    
4740
       !---- Local variables ----!
4741
       character(len=20), dimension(30) :: label
4742
       integer                          :: j,ic,n,na,nb,nc,nd,npos
4743
       integer                          :: nl,nl2,iv,ik
4744
       integer, dimension(1)            :: ivet
4745
       real(kind=cp)                    :: fac_0,fac_1
4746
       real(kind=cp),dimension(1)       :: vet
4747

    
4748
       call init_err_refcodes()
4749

    
4750
       nl=0
4751
       call getword(line,label,ic)
4752
       if (ic < 2) then
4753
          err_refcodes=.true.
4754
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
4755
          return
4756
       end if
4757

    
4758
       !---- Set Futher Information----!
4759
       !---- Na is the number of atom on List
4760
       !---- Nb is the number of keys (Rx,Ry,Rz,Ix,Iy,Iz,MagPh...)
4761
       !---- Fac0 is the multiplier
4762
       !---- Nl is the number of refinement parameter
4763

    
4764
       npos=index(label(1),"_")
4765
       if (npos ==0) then
4766
          err_refcodes=.true.
4767
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4768
          return
4769
       end if
4770

    
4771
       na=0
4772
       do j=1,FmAtom%Natoms
4773
          if (u_case(FmAtom%atom(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4774
             na=j
4775
             exit
4776
          end if
4777
       end do
4778
       if (na == 0) then
4779
          err_refcodes=.true.
4780
          ERR_RefCodes_Mess=" Atom label not found for "//trim(line)
4781
          return
4782
       end if
4783

    
4784
       nb=0
4785
       do j=1,mNcode
4786
          if (u_case(label(1)(1:npos))==u_case(trim(mcode_nam(j)))) then
4787
             nb=j
4788
             exit
4789
          end if
4790
       end do
4791
       if (nb == 0) then
4792
          err_refcodes=.true.
4793
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
4794
          return
4795
       end if
4796

    
4797
       !---- Get im, number of the magnetic matrices/irrep set
4798
       ik=FmAtom%atom(na)%nvk
4799
       
4800
       select case (nb)
4801
          case ( 1:3)
4802
             fac_0=FmAtom%atom(na)%mSkR(nb,ik)
4803
                nl=FmAtom%atom(na)%lSkR(nb,ik)
4804
          case ( 4:6)
4805
             fac_0=FmAtom%atom(na)%mSkI(nb-3,ik)
4806
                nl=FmAtom%atom(na)%lSkI(nb-3,ik)
4807
          case ( 7:9)
4808
             fac_0=FmAtom%atom(na)%mSkR(nb-6,ik)
4809
                nl=FmAtom%atom(na)%lSkR(nb-6,ik)
4810
          case (10:12)
4811
             fac_0=FmAtom%atom(na)%mSkI(nb-9,ik)
4812
                nl=FmAtom%atom(na)%lSkI(nb-9,ik)
4813
          case (13)
4814
             fac_0=FmAtom%atom(na)%mmphas(ik)
4815
                nl=FmAtom%atom(na)%lmphas(ik)
4816
           case (14:22)
4817
             fac_0=FmAtom%atom(na)%mbas(nb-13,ik)
4818
                nl=FmAtom%atom(na)%lbas(nb-13,ik)
4819
         case (23:)
4820
             err_refcodes=.true.
4821
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4822
             return
4823
       end select ! nb
4824

    
4825
       if (nb < mNcode) then
4826
          if (nl == 0) then
4827
             err_refcodes=.true.
4828
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
4829
             return
4830
          end if
4831
       end if
4832

    
4833
       !---- Set the rest elements in Contsraints ----!
4834
       n=1
4835
       do
4836
          n=n+1
4837
          if (n > ic) exit
4838

    
4839
          npos=index(label(n),"_")
4840
          if (npos ==0) then
4841
             err_refcodes=.true.
4842
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
4843
             return
4844
          end if
4845

    
4846
          nc=0
4847
          do j=1,FmAtom%Natoms
4848
             if (u_case(FmAtom%atom(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
4849
                nc=j
4850
                exit
4851
             end if
4852
          end do
4853
          if (nc == 0) then
4854
             err_refcodes=.true.
4855
             ERR_RefCodes_Mess=" Atom label not found for "//trim(label(n))
4856
             return
4857
          end if
4858

    
4859
          nd=0
4860
          do j=1,ncode
4861
             if (u_case(label(n)(1:npos))==u_case(trim(mcode_nam(j)))) then
4862
                nd=j
4863
                exit
4864
             end if
4865
          end do
4866
          if (nd == 0) then
4867
             err_refcodes=.true.
4868
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
4869
             return
4870
          end if
4871

    
4872
          !---- Is there a new multiplier?
4873
          n=n+1
4874
          call getnum(label(n),vet,ivet,iv)
4875
          if (iv == 1) then
4876
             fac_1=vet(1)
4877
          else
4878
             fac_1=fac_0
4879
             n=n-1
4880
          end if
4881

    
4882
          select case (nd)
4883
             case ( 1:3)
4884
                nl2=FmAtom%atom(nc)%lSkR(nd,ik)
4885
                call Delete_refCodes(nl2,FmAtom)
4886
                FmAtom%atom(nc)%mSkR(nd,ik)=fac_1
4887
                FmAtom%atom(nc)%lSkR(nd,ik)=nl
4888
             case ( 4:6)
4889
                nl2=FmAtom%atom(nc)%lSkI(nd-3,ik)
4890
                call Delete_refCodes(nl2,FmAtom)
4891
                FmAtom%atom(nc)%mSkI(nd-3,ik)=fac_1
4892
                FmAtom%atom(nc)%lSkI(nd-3,ik)=nl
4893
             case ( 7:9)
4894
                nl2=FmAtom%atom(nc)%lSkR(nd-6,ik)
4895
                call Delete_refCodes(nl2,FmAtom)
4896
                FmAtom%atom(nc)%mSkR(nd-6,ik)=fac_1
4897
                FmAtom%atom(nc)%lSkR(nd-6,ik)=nl
4898
             case (10:12)
4899
                nl2=FmAtom%atom(nc)%lSkI(nd-9,ik)
4900
                call Delete_refCodes(nl2,FmAtom)
4901
                FmAtom%atom(nc)%mSkI(nd-9,ik)=fac_1
4902
                FmAtom%atom(nc)%lSkI(nd-9,ik)=nl
4903
             case (13)
4904
                nl2=FmAtom%atom(nc)%lmphas(ik)
4905
                call Delete_refCodes(nl2,FmAtom)
4906
                FmAtom%atom(nc)%mmphas(ik)=fac_1
4907
                FmAtom%atom(nc)%lmphas(ik)=nl
4908
             case (14:22)
4909
                nl2=FmAtom%atom(nc)%lbas(nd-13,ik)
4910
                call Delete_refCodes(nl2,FmAtom)
4911
                FmAtom%atom(nc)%mbas(nd-13,ik)=fac_1
4912
                FmAtom%atom(nc)%lbas(nd-13,ik)=nl
4913
             case (23:)
4914
                err_refcodes=.true.
4915
                ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4916
                return
4917
          end select ! nb
4918

    
4919
          np_cons=np_cons+1
4920

    
4921
       end do
4922

    
4923
       return
4924
    End Subroutine Get_ConCodes_Line_FmAtom
4925

    
4926
    !!--++
4927
    !!--++ Subroutine Get_ConCodes_Line_Molcrys(Line,Molcrys)
4928
    !!--++    character(len=*),             intent(in)     :: Line
4929
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
4930
    !!--++
4931
    !!--++ Overloaded
4932
    !!--++ Get the Constraints relations
4933
    !!--++
4934
    !!--++ Update: March - 2005
4935
    !!
4936
    Subroutine Get_ConCodes_Line_Molcrys(Line,Molcrys)
4937
       !---- Arguments ----!
4938
       character(len=*),             intent(in)     :: Line
4939
       type(molecular_Crystal_type), intent(in out) :: MolCrys
4940

    
4941
       !---- Loval variables ----!
4942
       character(len=5)                 :: car
4943
       character(len=20), dimension(30) :: label
4944
       integer                          :: i,j,ic,n,na,naa,nb,nc,ncc,nd
4945
       integer                          :: npos, nposm, nmol1,nmol2
4946
       integer                          :: nl,nl2,iv
4947
       integer, dimension(1)            :: ivet
4948
       real(kind=cp)                    :: fac_0,fac_1
4949
       real(kind=cp),dimension(1)       :: vet
4950

    
4951
       call init_err_refcodes()
4952

    
4953
       call getword(line,label,ic)
4954
       if (ic < 2) then
4955
          err_refcodes=.true.
4956
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
4957
          return
4958
       end if
4959

    
4960
       !---- Set Father Information ----!
4961
       !---- Na is the number of atom on List
4962
       !---- Nb is the key (X,Y,Z,Occ,...)
4963
       !---- Fac0 is the multiplier
4964
       !---- Nl is the number of refinement parameter
4965
       npos=index(label(1),"_")
4966
       if (npos ==0) then
4967
          err_refcodes=.true.
4968
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4969
          return
4970
       end if
4971
       nposm=index(u_case(label(1)),"MOL")
4972
       if (nposm /= 0) then
4973
          car=adjustl(label(1)(nposm+3:))
4974
          if (car(1:2) == "  ") then
4975
             nmol1=0
4976
          else
4977
             read(unit=car,fmt="(i2)") nmol1
4978
          end if
4979
       else
4980
          nmol1=-1
4981
       end if
4982

    
4983
       na=0
4984
       do j=1,molcrys%n_free
4985
          if (u_case(molcrys%atm(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4986
             na=j
4987
             exit
4988
          end if
4989
       end do
4990
       if (na == 0) then
4991
          do i=1,molcrys%n_mol
4992
             do j=1,molcrys%mol(i)%natoms
4993
                if (u_case(molcrys%mol(i)%Atname(j)) == u_case(label(1)(npos+1:npos+6))) then
4994
                   if (j > 1) then
4995
                      na=molcrys%n_free+sum(molcrys%mol(1:j-1)%natoms)+j
4996
                   else
4997
                      na=molcrys%n_free+j
4998
                   end if
4999
                   exit
5000
                end if
5001
             end do
5002
          end do
5003
       end if
5004

    
5005
       nb=0
5006
       do j=1,ncode
5007
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
5008
             nb=j
5009
             exit
5010
          end if
5011
       end do
5012
       if (nb == 0) then
5013
          err_refcodes=.true.
5014
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
5015
          return
5016
       end if
5017

    
5018
       !---- Checking ----!
5019
       if (nb <= 12 .and. na==0) then
5020
          err_refcodes=.true.
5021
          ERR_RefCodes_Mess="Incompatible option for parameter name: "//trim(label(1))
5022
          return
5023
       end if
5024

    
5025
       nl=0
5026
       select case (nb)
5027
          case ( 1:3)
5028
             !---- X_, Y_, Z_ ----!
5029
             if (na <= molcrys%n_free) then
5030
                fac_0=molcrys%atm(na)%mx(nb)
5031
                   nl=molcrys%atm(na)%lx(nb)
5032
             else
5033
                naa=na-molcrys%n_free
5034
                do i=1,molcrys%n_mol
5035
                   if (naa > molcrys%mol(i)%natoms) then
5036
                      naa=naa-molcrys%mol(i)%natoms
5037
                      cycle
5038
                   end if
5039
                   fac_0=molcrys%mol(i)%mI_coor(nb,naa)
5040
                      nl=molcrys%mol(i)%lI_coor(nb,naa)
5041
                end do
5042
             end if
5043

    
5044
          case ( 4)
5045
             !---- Biso_ ----!
5046
             if (na <= molcrys%n_free) then
5047
                fac_0=molcrys%atm(na)%mbiso
5048
                   nl=molcrys%atm(na)%lbiso
5049
             else
5050
                naa=na-molcrys%n_free
5051
                do i=1,molcrys%n_mol
5052
                   if (naa > molcrys%mol(i)%natoms) then
5053
                      naa=naa-molcrys%mol(i)%natoms
5054
                      cycle
5055
                   end if
5056
                   fac_0=molcrys%mol(i)%mbiso(naa)
5057
                      nl=molcrys%mol(i)%lbiso(naa)
5058
                end do
5059
             end if
5060

    
5061
          case ( 5)
5062
             !---- Occ_ ----!
5063
             if (na <= molcrys%n_free) then
5064
                fac_0=molcrys%atm(na)%mocc
5065
                   nl=molcrys%atm(na)%locc
5066
             else
5067
                naa=na-molcrys%n_free
5068
                do i=1,molcrys%n_mol
5069
                   if (naa > molcrys%mol(i)%natoms) then
5070
                      naa=naa-molcrys%mol(i)%natoms
5071
                      cycle
5072
                   end if
5073
                   fac_0=molcrys%mol(i)%mocc(naa)
5074
                      nl=molcrys%mol(i)%locc(naa)
5075
                end do
5076
             end if
5077

    
5078
          case ( 6:11)
5079
             !---- B11_, ..., B23_ ----!
5080
             if (na <= molcrys%n_free) then
5081
                fac_0=molcrys%atm(na)%mu(nb-5)
5082
                   nl=molcrys%atm(na)%lu(nb-5)
5083
             else
5084
                err_refcodes=.true.
5085
                ERR_RefCodes_Mess="Option no valid"
5086
                return
5087
             end if
5088

    
5089
          case (12)
5090
             !---- Banis_ ----!
5091
             if (na <= molcrys%n_free) then
5092
                fac_0=molcrys%atm(na)%mu(1)
5093
             else
5094
                err_refcodes=.true.
5095
                ERR_RefCodes_Mess="Option no valid"
5096
                return
5097
             end if
5098

    
5099
          case (13:15)
5100
             !---- Xc_, Yc_, Zc_ ----!
5101
             select case (nmol1)
5102
                case (-1)
5103
                   err_refcodes=.true.
5104
                   ERR_RefCodes_Mess="Option no valid"
5105
                   return
5106

    
5107
                case (0)
5108
                   fac_0=molcrys%mol(1)%mxcentre(nb-12)
5109
                   nl=molcrys%mol(1)%lxcentre(nb-12)
5110

    
5111
                case (1:)
5112
                   fac_0=molcrys%mol(nmol1)%mxcentre(nb-12)
5113
                   nl=molcrys%mol(nmol1)%lxcentre(nb-12)
5114
             end select
5115

    
5116
          case (16:18)
5117
             !---- Theta_ , Phi_, Chi_ ----!
5118
             select case (nmol1)
5119
                case (-1)
5120
                   err_refcodes=.true.
5121
                   ERR_RefCodes_Mess="Option no valid"
5122
                   return
5123

    
5124
                case (0)
5125
                   fac_0=molcrys%mol(1)%mOrient(nb-15)
5126
                   nl=molcrys%mol(1)%lOrient(nb-15)
5127

    
5128
                case (1)
5129
                   fac_0=molcrys%mol(nmol1)%mOrient(nb-15)
5130
                   nl=molcrys%mol(nmol1)%lOrient(nb-15)
5131
             end select
5132

    
5133
          case (19:21)
5134
             !!! Not yet Implemented !!!
5135
       end select ! nb
5136

    
5137
       if (nb < ncode) then
5138
          if (nl == 0) then
5139
             err_refcodes=.true.
5140
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
5141
             return
5142
          end if
5143
       end if
5144

    
5145
       !---- Set the rest elements in Constraints ----!
5146
       n=1
5147
       do
5148
          n=n+1
5149
          if (n > ic) exit
5150

    
5151
          npos=index(label(n),"_")
5152
          if (npos ==0) then
5153
             err_refcodes=.true.
5154
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
5155
             return
5156
          end if
5157
          nposm=index(u_case(label(n)),"MOL")
5158
          if (nposm /= 0) then
5159
             car=adjustl(label(n)(nposm+3:))
5160
             if (car(1:2) == "  ") then
5161
                nmol2=0
5162
             else
5163
                read(unit=car,fmt="(i2)") nmol2
5164
             end if
5165
          else
5166
             nmol2=-1
5167
          end if
5168

    
5169
          nc=0
5170
          do j=1,molcrys%n_free
5171
             if (u_case(molcrys%atm(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
5172
                nc=j
5173
                exit
5174
             end if
5175
          end do
5176
          if (nc == 0) then
5177
             do i=1,molcrys%n_mol
5178
                do j=1,molcrys%mol(i)%natoms
5179
                   if (u_case(molcrys%mol(i)%Atname(j)) == u_case(label(n)(npos+1:npos+6))) then
5180
                      if (j > 1) then
5181
                         nc=molcrys%n_free+sum(molcrys%mol(1:j-1)%natoms)+j
5182
                      else
5183
                         nc=molcrys%n_free+j
5184
                      end if
5185
                      exit
5186
                   end if
5187
                end do
5188
             end do
5189
          end if
5190

    
5191
          nd=0
5192
          do j=1,ncode
5193
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
5194
                nd=j
5195
                exit
5196
             end if
5197
          end do
5198
          if (nd == 0) then
5199
             err_refcodes=.true.
5200
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
5201
             return
5202
          end if
5203

    
5204
          !---- Checking ----!
5205
          if (nd <= 12 .and. nc==0) then
5206
             err_refcodes=.true.
5207
             ERR_RefCodes_Mess="Incompatible option for parameter name: "//trim(label(n))
5208
             return
5209
          end if
5210

    
5211
          !---- Is there a new multiplier? ----!
5212
          n=n+1
5213
          call getnum(label(n),vet,ivet,iv)
5214
          if (iv == 1) then
5215
             fac_1=vet(1)
5216
          else
5217
             fac_1=fac_0
5218
             n=n-1
5219
          end if
5220

    
5221
          select case (nd)
5222
             case ( 1:3)
5223
                !---- X_, Y_, Z_ ----!
5224
                if (nc <= molcrys%n_free) then
5225
                   nl2=molcrys%atm(nc)%lx(nd)
5226
                   call Delete_RefCodes(nl2,molcrys)
5227
                   molcrys%atm(nc)%mx(nd)=fac_1
5228
                   molcrys%atm(nc)%lx(nd)=nl
5229
                else
5230
                   ncc=nc-molcrys%n_free
5231
                   do i=1,molcrys%n_mol
5232
                      if (ncc > molcrys%mol(i)%natoms) then
5233
                         ncc=ncc-molcrys%mol(i)%natoms
5234
                         cycle
5235
                      end if
5236
                      nl2=molcrys%mol(i)%lI_coor(nd,ncc)
5237
                      call Delete_RefCodes(nl2,molcrys)
5238
                      molcrys%mol(i)%mI_coor(nd,ncc)=fac_1
5239
                      molcrys%mol(i)%lI_coor(nd,ncc)=nl
5240
                   end do
5241
                end if
5242

    
5243
             case ( 4)
5244
                !---- Biso_ ----!
5245
                if (nc <= molcrys%n_free) then
5246
                   nl2=molcrys%atm(nc)%lbiso
5247
                   call Delete_RefCodes(nl2,molcrys)
5248
                   molcrys%atm(nc)%mbiso=fac_1
5249
                   molcrys%atm(nc)%lbiso=nl
5250
                else
5251
                   ncc=nc-molcrys%n_free
5252
                   do i=1,molcrys%n_mol
5253
                      if (ncc > molcrys%mol(i)%natoms) then
5254
                         ncc=ncc-molcrys%mol(i)%natoms
5255
                         cycle
5256
                      end if
5257
                      nl2=molcrys%mol(i)%lbiso(ncc)
5258
                      call Delete_RefCodes(nl2,molcrys)
5259
                      molcrys%mol(i)%mbiso(ncc)=fac_1
5260
                      molcrys%mol(i)%lbiso(ncc)=nl
5261
                   end do
5262
                end if
5263

    
5264
             case ( 5)
5265
                !---- Occ_ ----!
5266
                if (nc <= molcrys%n_free) then
5267
                   nl2=molcrys%atm(nc)%locc
5268
                   call Delete_RefCodes(nl2,molcrys)
5269
                   molcrys%atm(nc)%mocc=fac_1
5270
                   molcrys%atm(nc)%locc=nl
5271
                else
5272
                   ncc=nc-molcrys%n_free
5273
                   do i=1,molcrys%n_mol
5274
                      if (ncc > molcrys%mol(i)%natoms) then
5275
                         ncc=ncc-molcrys%mol(i)%natoms
5276
                         cycle
5277
                      end if
5278
                      nl2=molcrys%mol(i)%locc(ncc)
5279
                      call Delete_RefCodes(nl2,molcrys)
5280
                      molcrys%mol(i)%mocc(ncc)=fac_1
5281
                      molcrys%mol(i)%locc(ncc)=nl
5282
                   end do
5283
                end if
5284

    
5285
             case ( 6:11)
5286
                !---- B11_, ...., B23_ ----!
5287
                if (nc <= molcrys%n_free) then
5288
                   nl2=molcrys%atm(nc)%lu(nd-5)
5289
                   call Delete_RefCodes(nl2,molcrys)
5290
                   molcrys%atm(nc)%mu(nd-5)=fac_1
5291
                   molcrys%atm(nc)%lu(nd-5)=nl
5292
                else
5293
                   err_refcodes=.true.
5294
                   ERR_RefCodes_Mess="Option no valid"
5295
                   return
5296
                end if
5297

    
5298
             case (12)
5299
                !---- Banis_ ----!
5300
                if (nc <= molcrys%n_free) then
5301
                   do j=1,6
5302
                      nl2=molcrys%atm(nc)%lu(j)
5303
                      call Delete_RefCodes(nl2,molcrys)
5304
                      molcrys%atm(nc)%mu(j)=fac_1
5305
                      molcrys%atm(nc)%lu(j)=molcrys%atm(na)%lu(j)
5306
                   end do
5307
                   np_cons=np_cons+5
5308
                else
5309
                   err_refcodes=.true.
5310
                   ERR_RefCodes_Mess="Option no valid"
5311
                   return
5312
                end if
5313

    
5314
             case (13:15)
5315
                select case (nmol2)
5316
                   case (-1)
5317
                      err_refcodes=.true.
5318
                      ERR_RefCodes_Mess="Option no valid"
5319
                      return
5320

    
5321
                   case (0)
5322
                      do i=1,molcrys%n_mol
5323
                         nl2=molcrys%mol(i)%lxcentre(nc-12)
5324
                         call Delete_RefCodes(nl2,molcrys)
5325
                         molcrys%mol(i)%mxcentre(nc-12)=fac_1
5326
                         molcrys%mol(i)%lxcentre(nc-12)=nl
5327
                      end do
5328
                      np_cons=np_cons+(i-1)
5329

    
5330
                   case (1:)
5331
                      nl2=molcrys%mol(nmol2)%lxcentre(nc-12)
5332
                      call Delete_RefCodes(nl2,molcrys)
5333
                      molcrys%mol(nmol2)%mxcentre(nc-12)=fac_1
5334
                      molcrys%mol(nmol2)%lxcentre(nc-12)=nl
5335
                end select
5336

    
5337
             case (16:18)
5338
                select case (nmol2)
5339
                   case (-1)
5340
                      err_refcodes=.true.
5341
                      ERR_RefCodes_Mess="Option no valid"
5342
                      return
5343

    
5344
                   case (0)
5345
                      do i=1,molcrys%n_mol
5346
                         nl2=molcrys%mol(i)%lorient(nc-15)
5347
                         call Delete_RefCodes(nl2,molcrys)
5348
                         molcrys%mol(i)%morient(nc-15)=fac_1
5349
                         molcrys%mol(i)%lorient(nc-15)=nl
5350
                      end do
5351
                      np_cons=np_cons+(i-1)
5352

    
5353
                   case (1:)
5354
                      nl2=molcrys%mol(nmol2)%lorient(nc-15)
5355
                      call Delete_RefCodes(nl2,molcrys)
5356
                      molcrys%mol(nmol2)%morient(nc-15)=fac_1
5357
                      molcrys%mol(nmol2)%lorient(nc-15)=nl
5358
                end select
5359

    
5360
             case (19:21)
5361
                !!! Not yet implemented !!!
5362

    
5363
          end select ! nb
5364
          np_cons=np_cons+1
5365

    
5366
       end do
5367

    
5368
       return
5369
    End Subroutine Get_ConCodes_Line_Molcrys
5370

    
5371
    !!--++
5372
    !!--++ Subroutine Get_ConCodes_Line_Molec(Line,Molec)
5373
    !!--++    character(len=*),    intent(in)     :: Line
5374
    !!--++    type(molecule_type), intent(in out) :: Molec
5375
    !!--++
5376
    !!--++ Overloaded
5377
    !!--++ Get the Constraints relations
5378
    !!--++
5379
    !!--++ Update: March - 2005
5380
    !!
5381
    Subroutine Get_ConCodes_Line_Molec(Line,Molec)
5382
       !---- Arguments ----!
5383
       character(len=*),    intent(in)     :: Line
5384
       type(molecule_type), intent(in out) :: Molec
5385

    
5386
       !---- Loval variables ----!
5387
       character(len=20), dimension(30) :: label
5388
       integer                          :: j,ic,na,nb,nc,nd,npos!,i,naa,ncc
5389
       integer                          :: n,nl,nl2,iv
5390
       integer, dimension(1)            :: ivet
5391
       real(kind=cp)                    :: fac_0,fac_1
5392
       real(kind=cp),dimension(1)       :: vet
5393

    
5394
       call init_err_refcodes()
5395

    
5396
       call getword(line,label,ic)
5397
       if (ic < 2) then
5398
          err_refcodes=.true.
5399
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
5400
          return
5401
       end if
5402

    
5403
       !---- Set Father Information ----!
5404
       !---- Na is the number of atom on List
5405
       !---- Nb is the key (X,Y,Z,Occ,...)
5406
       !---- Fac0 is the multiplier
5407
       !---- Nl is the number of refinement parameter
5408
       npos=index(label(1),"_")
5409
       if (npos ==0) then
5410
          err_refcodes=.true.
5411
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
5412
          return
5413
       end if
5414

    
5415
       na=0
5416
       do j=1,molec%natoms
5417
          if (u_case(molec%Atname(j)) == u_case(label(1)(npos+1:npos+6))) then
5418
             na=j
5419
             exit
5420
          end if
5421
       end do
5422

    
5423
       nb=0
5424
       do j=1,ncode
5425
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
5426
             nb=j
5427
             exit
5428
          end if
5429
       end do
5430
       if (nb == 0) then
5431
          err_refcodes=.true.
5432
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
5433
          return
5434
       end if
5435

    
5436
       !---- Checking ----!
5437
       if (nb < 6 .and. na == 0) then
5438
          err_refcodes=.true.
5439
          ERR_RefCodes_Mess="Incompatible relation: "//trim(label(1))
5440
          return
5441
       end if
5442

    
5443
       nl=0
5444
       select case (nb)
5445
          case ( 1:3)
5446
             !---- X_, Y_, Z_ ----!
5447
             fac_0=molec%mI_coor(nb,na)
5448
                nl=molec%lI_coor(nb,na)
5449
          case ( 4)
5450
             !---- Biso_ ----!
5451
             fac_0=molec%mbiso(na)
5452
                nl=molec%lbiso(na)
5453
          case ( 5)
5454
             !---- Occ_ ----!
5455
             fac_0=molec%mocc(na)
5456
                nl=molec%locc(na)
5457
          case ( 6:12)
5458
             !---- Anisotropic Parameters ----!
5459
             err_refcodes=.true.
5460
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
5461
             return
5462
          case (13:15)
5463
             !---- Xc_, Yc_, Zc_ ----!
5464
             fac_0=molec%mxcentre(nb-12)
5465
                nl=molec%lxcentre(nb-12)
5466
          case (16:18)
5467
             !---- Theta_, Phi_, Chi_ ----!
5468
             fac_0=molec%mxcentre(nb-15)
5469
                nl=molec%lxcentre(nb-15)
5470
          case (19:21)
5471
             !!! Not Yet Implemented !!!
5472
       end select ! nb
5473

    
5474
       if (nb < ncode) then
5475
          if (nl == 0) then
5476
             err_refcodes=.true.
5477
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
5478
             return
5479
          end if
5480
       end if
5481

    
5482
       !---- Set Others ----!
5483
       n=1
5484
       do
5485
          n=n+1
5486
          if (n > ic) exit
5487

    
5488
          npos=index(label(n),"_")
5489
          if (npos ==0) then
5490
             err_refcodes=.true.
5491
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
5492
             return
5493
          end if
5494

    
5495
          nc=0
5496
          do j=1,molec%natoms
5497
             if (u_case(molec%Atname(j)) == u_case(label(n)(npos+1:npos+6))) then
5498
                nc=j
5499
                exit
5500
             end if
5501
          end do
5502

    
5503
          nd=0
5504
          do j=1,ncode
5505
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
5506