Project

General

Profile

CFML_Refcodes.f90

Zaharko, Oksana, 22/01/2012 10:29 AM

Download (342 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
!!----
12
!!----
13
!!---- DEPENDENCIES
14
!!----
15
!!----
16
!!---- VARIABLES
17
!!--..    Types
18
!!----    ANGLE_RESTRAINT_TYPE
19
!!----    DISTANCE_RESTRAINT_TYPE
20
!!----    TORSION_RESTRAINT_TYPE
21
!!--..
22
!!----    ANG_REST
23
!!----    CODE_NAM
24
!!----    mCODE_NAM                    !NEW 12/11
25
!!----    DIS_REST
26
!!----    ERR_REFCODES
27
!!----    ERR_REFCODES_MESS
28
!!--++    NCODE                        [Private]
29
!!--++    NKEY                         [Private]
30
!!--++    mNCODE                       [Private] !NEW 12/11
31
!!--++    mNKEY                        [Private] !NEW 12/11
32
!!----    KEY_CODE
33
!!----    KEY_mCODE                    !NEW 12/11
34
!!----    NP_CONS
35
!!----    NP_MAX
36
!!----    NP_REFI
37
!!----    NP_REST_ANG
38
!!----    NP_REST_DIS
39
!!----    NP_REST_TOR
40
!!----    TORSION_REST
41
!!----    V_BCON
42
!!----    V_BOUNDS
43
!!----    V_LIST
44
!!----    V_NAME
45
!!----    V_VEC
46
!!----    V_SHIFT
47
!!----
48
!!---- PUBLIC PROCEDURES
49
!!----    Functions:
50
!!----
51
!!----    Subroutines:
52
!!----       ALLOCATE_RESTPARAM
53
!!----       ALLOCATE_VPARAM
54
!!--++       DELETE_REFCODES             [Private]
55
!!--++       DELETE_REFCODES_FATOM       [Overloaded]
56
!!--++       DELETE_REFCODES_MOLCRYS     [Overloaded]
57
!!--++       DELETE_REFCODES_MOLEC       [Overloaded]
58
!!--++       
59
!!--++       FILL_REFCODES_FATOM         [Overloaded]
60
!!--++       FILL_REFCODES_MOLCRYS       [Overloaded]
61
!!--++       FILL_REFCODES_MOLEC         [Overloaded]
62
!!--++       GET_ATOMBET_CTR             [Private]
63
!!--++       GET_ATOMPOS_CTR             [Private]
64
!!--++       GET_CONCODES_LINE           [Private]
65
!!--++       GET_CONCODES_LINE_FATOM     [Overloaded]
66
!!--++       GET_CONCODES_LINE_MOLCRYS   [Overloaded]
67
!!--++       GET_CONCODES_LINE_MOLEC     [Overloaded]
68
!!--++       GET_REFCODES_LINE           [Private]
69
!!--++       GET_REFCODES_LINE_FATOM     [Overloaded]
70
!!--++       GET_REFCODES_LINE_MOLCRYS   [Overloaded]
71
!!--++       GET_REFCODES_LINE_MOLEC     [Overloaded]
72
!!----       GET_RESTANG_LINE
73
!!----       GET_RESTDIS_LINE
74
!!----       GET_RESTTOR_LINE
75
!!----       INIT_ERR_REFCODES
76
!!----       INIT_REFCODES
77
!!--++       INIT_REFCODES_FATOM         [Overloaded]
78
!!--++       INIT_REFCODES_FmATOM        [Overloaded] !NEW 12/11
79
!!--++       INIT_REFCODES_MOLCRYS       [Overloaded]
80
!!--++       INIT_REFCODES_MOLEC         [Overloaded]
81
!!----       READ_REFCODES_FILE
82
!!--++       READ_REFCODES_FILE_FATOM    [Overloaded]
83
!!--++       READ_REFCODES_FILE_FmATOM   [Overloaded] !NEW 12/11
84
!!--++       READ_REFCODES_FILE_MOLCRYS  [Overloaded]
85
!!--++       READ_REFCODES_FILE_MOLEC    [Overloaded]
86
!!--++       SPLIT_OPERATIONS            [Private]
87
!!--++       SPLIT_mOPERATIONS           [Private]    !NEW 12/11
88
!!----       VSTATE_TO_ATOMSPAR
89
!!--++       VSTATE_TO_ATOMSPAR_FATOM    [Overloaded]
90
!!--++       VSTATE_TO_ATOMSPAR_FmATOM   [Overloaded] !NEW 12/11
91
!!--++       VSTATE_TO_ATOMSPAR_MOLCRYS  [Overloaded]
92
!!--++       VSTATE_TO_ATOMSPAR_MOLEC    [Overloaded]
93
!!----       WRITE_INFO_REFCODES
94
!!--++       WRITE_INFO_REFCODES_FATOM   [Overloaded]
95
!!--++       WRITE_INFO_REFCODES_FmATOM  [Overloaded] !NEW 12/11
96
!!--++       WRITE_INFO_REFCODES_MOLCRYS [Overloaded]
97
!!--++       WRITE_INFO_REFCODES_MOLEC   [Overloaded]
98
!!----       WRITE_INFO_REFPARAMS
99
!!----       WRITE_RESTRAINTS_OBSCALC
100
!!----
101
!!
102
 Module CFML_Keywords_Code_Parser
103
    !---- Modules ----!
104
    Use CFML_GlobalDeps,                only: cp
105
    Use CFML_Math_General,              only: Sort
106
    Use CFML_String_Utilities,          only: Cutst, U_Case, L_Case, Getword, GetNum
107
    Use CFML_Crystallographic_Symmetry, only: Space_Group_Type, Get_Stabilizer, Symmetry_Symbol,   &
108
                                              Sym_B_Relations, Read_SymTrans_Code, Get_SymSymb,Read_Xsym
109
    Use CFML_Atom_TypeDef,              only: Atom_list_Type, mAtom_list_Type
110
    Use CFML_Molecular_Crystals,        only: Molecule_Type, Molecular_Crystal_Type
111
    Use CFML_IO_Formats,                only: File_List_Type
112
    Use CFML_Magnetic_Symmetry
113
    Use CFML_Math_3D,                   only: Get_Spheric_Coord
114

    
115
    !---- Variables ----!
116
    implicit none
117

    
118
    private
119

    
120
    !---- List of public functions ----!
121

    
122
    !---- List of public subroutines ----!
123
    public :: Allocate_VParam, Init_RefCodes, Read_RefCodes_File, VState_to_AtomsPar,  &
124
              Write_Info_RefCodes, Get_RestAng_Line, Get_RestDis_Line, Get_RestTor_Line, &
125
              Allocate_RestParam, Write_Restraints_ObsCalc, Init_Err_RefCodes, &
126
              Write_Info_RefParams
127

    
128
    !---- List of private functions ----!
129

    
130
    !---- List of private subroutines ----!
131
    private :: Delete_RefCodes,   &
132
               Delete_RefCodes_FAtom, Delete_RefCodes_FmAtom, Delete_RefCodes_Molcrys,          &
133
               Delete_RefCodes_Molec, &
134
               Fill_RefCodes,     &
135
               Fill_RefCodes_FAtom, Fill_RefCodes_FmAtom, Fill_RefCodes_Molcrys,                &
136
               Fill_RefCodes_Molec,   &
137
               Get_AtomBet_Ctr, Get_Atompos_Ctr,                                                &
138
               Get_ConCodes_Line, &
139
               Get_ConCodes_Line_FAtom, Get_ConCodes_Line_FmAtom, Get_ConCodes_Line_Molcrys,    &
140
               Get_ConCodes_Line_Molec, Get_RefCodes_Line, &
141
               Get_RefCodes_Line_FAtom, Get_RefCodes_Line_FmAtom, Get_RefCodes_Line_Molcrys,    &
142
               Get_RefCodes_Line_Molec, &
143
               Init_RefCodes_FAtom, Init_RefCodes_FmAtom, Init_RefCodes_Molcrys,                &
144
               Init_RefCodes_Molec, &
145
               Read_RefCodes_File_FAtom, Read_RefCodes_File_FmAtom, Read_RefCodes_File_Molcrys, &
146
               Read_RefCodes_File_Molec, Split_Operations, Split_mOperations, &
147
               VState_to_AtomsPar_FAtom, VState_to_AtomsPar_FmAtom, VState_to_AtomsPar_Molcrys, &
148
               VState_to_AtomsPar_Molec, &
149
               Write_Info_RefCodes_FAtom,Write_Info_RefCodes_FmAtom,Write_Info_RefCodes_Molcrys,& 
150
               Write_Info_RefCodes_Molec
151

    
152
    !---- Definitions ----!
153

    
154
    !!----
155
    !!---- TYPE :: ANGLE_RESTRAINT_TYPE
156
    !!--..
157
    !!---- Type, public :: Angle_Restraint_Type
158
    !!----    real(kind=cp)                 :: AObs
159
    !!----    real(kind=cp)                 :: ACalc
160
    !!----    real(kind=cp)                 :: Sigma
161
    !!----    integer,dimension(3)          :: P
162
    !!----    character(len=8),dimension(2) :: STCode
163
    !!---- End Type Angle_Restraint_Type
164
    !!----
165
    !!---- Update: April - 2005
166
    !!
167
    Type, public :: Angle_Restraint_Type
168
      real(kind=cp)                 :: AObs
169
      real(kind=cp)                 :: ACalc
170
      real(kind=cp)                 :: Sigma
171
      integer,dimension(3)          :: P
172
      character(len=8),dimension(2) :: STCode
173
    End Type Angle_Restraint_Type
174

    
175
    !!----
176
    !!---- TYPE :: DISTANCE_RESTRAINT_TYPE
177
    !!--..
178
    !!---- Type, public :: Distance_Restraint_Type
179
    !!----    real(kind=cp)        :: DObs
180
    !!----    real(kind=cp)        :: DCalc
181
    !!----    real(kind=cp)        :: Sigma
182
    !!----    integer,dimension(2) :: P
183
    !!----    character(len=8)     :: STCode    ! _N.ABC
184
    !!---- End Type Distance_Restraint_Type
185
    !!----
186
    !!---- Update: April - 2005
187
    !!
188
    Type, public :: Distance_Restraint_Type
189
      real(kind=cp)        :: DObs
190
      real(kind=cp)        :: DCalc
191
      real(kind=cp)        :: Sigma
192
      integer,dimension(2) :: P
193
      character(len=8)     :: STCode
194
    End Type Distance_Restraint_Type
195

    
196
    !!----
197
    !!---- TYPE :: TORSION_RESTRAINT_TYPE
198
    !!--..
199
    !!---- Type, public :: Torsion_Restraint_Type
200
    !!----    real(kind=cp)                 :: TObs
201
    !!----    real(kind=cp)                 :: TCalc
202
    !!----    real(kind=cp)                 :: Sigma
203
    !!----    integer,dimension(4)          :: P
204
    !!----    character(len=8),dimension(3) :: STCode
205
    !!---- End Type Torsion_Restraint_Type
206
    !!----
207
    !!---- Update: April - 2005
208
    !!
209
    Type, public :: Torsion_Restraint_Type
210
      real(kind=cp)                 :: TObs
211
      real(kind=cp)                 :: TCalc
212
      real(kind=cp)                 :: Sigma
213
      integer,dimension(4)          :: P
214
      character(len=8),dimension(3) :: STCode
215
    End Type Torsion_Restraint_Type
216

    
217
    !!----
218
    !!---- ANG_REST
219
    !!----    type(Angle_Restraint_Type), public, dimension(:), allocatable :: Ang_rest
220
    !!----
221
    !!---- Relations for Angle Restraints
222
    !!----
223
    !!---- Update: March - 2005
224
    !!
225
    type(Angle_Restraint_Type), public, dimension(:), allocatable :: Ang_Rest
226

    
227
    !!--++
228
    !!--++ NCODE
229
    !!--++    integer, private :: NCode
230
    !!--++
231
    !!--++    Number of Code variables
232
    !!--++
233
    !!--++ Update: March - 2005
234
    !!
235
    integer, private, parameter :: NCode=21
236

    
237
    !!----
238
    !!---- CODE_NAM
239
    !!----    character(len=*), dimension(NCode), public, parameter :: Code_Nam
240
    !!----
241
    !!----    Variable for treatement Codes
242
    !!--..    21 debe ser igual a NCode
243
    !!----
244
    !!---- Update: March - 2005
245
    !!
246
    character(len=*), dimension(NCode), public, parameter :: Code_Nam=(/ "X_    ","Y_    ","Z_    ",       &
247
                                                                         "B_    ","Occ_  ","B11_  ",       &
248
                                                                         "B22_  ","B33_  ","B12_  ",       &
249
                                                                         "B13_  ","B23_  ","Bns_  ",       &
250
                                                                         "Xc_   ","Yc_   ","Zc_   ",       &
251
                                                                         "Theta_","Phi_  ","Chi_  ",       &
252
                                                                         "Th_L_ ","Th_T_ ","Th_S_ "/)
253
    integer, private, parameter :: mNCode=25
254

    
255
    !!----
256
    !!---- mCODE_NAM
257
    !!----    character(len=*), dimension(mNCode), public, parameter :: mCode_Nam
258
    !!----
259
    !!----    Variable for treatement Codes
260
    !!----    mag clone of CODE_NAM
261
    !!---- Created: December - 2011
262
    !!
263
    character(len=*),dimension(mNCode), public, parameter :: mCode_Nam=(/"Rx_   ","Ry_   ","Rz_   ",&
264
                                                                         "Ix_   ","Iy_   ","Iz_   ",&
265
                                                                         "Rm_   ","Rphi_ ","Rth_  ",&
266
                                                                         "Im_   ","Iphi_ ","Ith_  ",&
267
                                                                         "MagPh_",               &
268
                                                                         "C1_   ","C2_   ","C3_   ",   &
269
                                                                         "C4_   ","C5_   ","C6_   ",   &
270
                                                                         "C7_   ","C8_   ","C9_   ",   &
271
                                                                         "C10_  ","C11_  ","C12_  "/)
272

    
273
    !!----
274
    !!---- DIS_REST
275
    !!----    type(Distance_Restraint_Type), public, dimension(:), allocatable :: Dis_Rest
276
    !!----
277
    !!---- Relations for Angle Restraints
278
    !!----
279
    !!---- Update: March - 2005
280
    !!
281
    type(Distance_Restraint_Type), public, dimension(:), allocatable :: Dis_Rest
282

    
283
    !!----
284
    !!---- ERR_REFCODES
285
    !!----    logical, public :: Err_RefCodes
286
    !!----
287
    !!----    Error variable
288
    !!----
289
    !!---- Update: March - 2005
290
    !!
291
    logical, public :: Err_RefCodes = .false.
292

    
293
    !!----
294
    !!---- ERR_REFCODES_MESS
295
    !!----    character(len=150), public :: ERR_RefCodes_Mess
296
    !!----
297
    !!----    Error variable messages
298
    !!----
299
    !!---- Update: March - 2005
300
    !!
301
    character(len=150), public :: ERR_RefCodes_Mess = " "
302

    
303
    !!--++
304
    !!--++ NKEY
305
    !!--++    integer, public :: NKey
306
    !!--++
307
    !!--++    Number of Keys variables
308
    !!--++
309
    !!--++ Update: March - 2005
310
    !!
311
    integer, private, parameter :: NKey=8
312

    
313
    !!----
314
    !!---- KEY_CODE
315
    !!----    character(len=*), dimension(NKey), public, parameter :: key_Code
316
    !!----
317
    !!----     Key codes defined in the module
318
    !!----
319
    !!---- Update: March - 2005
320
    !!
321
    character(len=*), dimension(Nkey), public, parameter :: Key_Code=(/ "XYZ ","OCC ","BIS ","BAN ", &
322
                                                                        "ALL ","CEN ","ORI ","THE "/)
323
    integer, private, parameter :: mNKey=3
324

    
325
    !!----
326
    !!---- KEY_mCODE
327
    !!----    character(len=*), dimension(mNKey), public, parameter :: key_mCode
328
    !!----
329
    !!----    Key codes defined in the module
330
    !!----    mag clone of KEY_CODE
331
    !!---- Created: December - 2011
332
    !!
333
    character(len=*), dimension(mNkey), public, parameter :: Key_mCode=(/"Rxyz","Ixyz","Mxyz"/)
334
    !!----
335
    !!---- NP_CONS
336
    !!----    integer, public :: NP_Cons
337
    !!----
338
    !!----    Number of Constraints relations
339
    !!----
340
    !!---- Update: March - 2005
341
    !!
342
    integer, public :: NP_Cons
343

    
344
    !!----
345
    !!---- NP_MAX
346
    !!----    integer, public :: NP_Max
347
    !!----
348
    !!----    Number of Maximum Parameters to Refine
349
    !!----
350
    !!---- Update: March - 2005
351
    !!
352
    integer, public :: NP_Max
353

    
354
    !!----
355
    !!---- NP_REFI
356
    !!----    integer, public :: NP_Refi
357
    !!----
358
    !!----    Number of Refinable Parameters
359
    !!----
360
    !!---- Update: March - 2005
361
    !!
362
    integer, public :: NP_Refi
363

    
364
    !!----
365
    !!---- NP_REST_ANG
366
    !!----    integer, public :: NP_Rest_Ang
367
    !!----
368
    !!----    Number of Angle Restraints relations
369
    !!----
370
    !!---- Update: March - 2005
371
    !!
372
    integer, public :: NP_Rest_Ang
373

    
374
    !!----
375
    !!---- NP_REST_DIS
376
    !!----    integer, public :: NP_Rest_Dis
377
    !!----
378
    !!----    Number of Distance Restraints relations
379
    !!----
380
    !!---- Update: March - 2005
381
    !!
382
    integer, public :: NP_Rest_Dis
383

    
384
    !!----
385
    !!---- NP_REST_TOR
386
    !!----    integer, public :: NP_Rest_Tor
387
    !!----
388
    !!----    Number of Torsion Restraints relations
389
    !!----
390
    !!---- Update: March - 2005
391
    !!
392
    integer, public :: NP_Rest_Tor
393

    
394
    !!----
395
    !!---- TOR_REST
396
    !!----    type(Torsion_Restraint_Type), public, dimension(:), allocatable :: Tor_Rest
397
    !!----
398
    !!---- Relations for Torsion Angle Restraints
399
    !!----
400
    !!---- Update: March - 2005
401
    !!
402
    type(Torsion_Restraint_Type), public, dimension(:), allocatable :: Tor_Rest
403

    
404
    !!----
405
    !!---- V_BCON
406
    !!----    integer, public, dimension(:), allocatable :: V_BCon
407
    !!----
408
    !!----    Vector of Boundary Conditions
409
    !!----
410
    !!---- Update: March - 2005
411
    !!
412
    integer, public, dimension(:), allocatable :: V_BCon
413

    
414
    !!----
415
    !!---- V_BOUNDS
416
    !!----    real(kind=cp), public, dimension(:,:), allocatable :: V_Bounds
417
    !!----
418
    !!----    Vector of Lower, Upper limits and Step for Parameters
419
    !!----
420
    !!---- Update: March - 2005
421
    !!
422
    real(kind=cp), public, dimension(:,:),  allocatable :: V_Bounds
423

    
424
    !!----
425
    !!---- V_LIST
426
    !!----    integer, public, dimension(:), allocatable :: V_List
427
    !!----
428
    !!----    Vector of Index point the atom order
429
    !!----
430
    !!---- Update: March - 2005
431
    !!
432
    integer, public, dimension(:),  allocatable :: V_List
433

    
434
    !!----
435
    !!---- V_NAME
436
    !!----    character(len=20), public, dimension(:), allocatable :: V_Name
437
    !!----
438
    !!----    Vector of  Name of Refinable Parameters
439
    !!----
440
    !!---- Update: March - 2005
441
    !!
442
    character(len=20), public, dimension(:), allocatable :: V_Name
443

    
444
    !!----
445
    !!---- V_VEC
446
    !!----    real(kind=cp), public, dimension(:), allocatable :: V_Vec
447
    !!----
448
    !!----    Vector of  Parameters
449
    !!----
450
    !!---- Update: March - 2005
451
    !!
452
    real(kind=cp), public, dimension(:),    allocatable :: V_Vec
453

    
454
    !!----
455
    !!---- V_SHIFT
456
    !!----    real(kind=cp), public, dimension(:), allocatable :: V_Shift
457
    !!----
458
    !!----    Vector of holding the shift of parameters
459
    !!----
460
    !!---- Update: March - 2005
461
    !!
462
    real(kind=cp), public, dimension(:),    allocatable :: V_Shift
463

    
464
    !---- Interfaces - Overloaded ----!
465
    Interface Delete_RefCodes
466
       Module Procedure Delete_RefCodes_FAtom
467
       Module Procedure Delete_RefCodes_FmAtom
468
       Module Procedure Delete_RefCodes_Molcrys
469
       Module Procedure Delete_RefCodes_Molec
470
    End Interface
471

    
472
    Interface Fill_RefCodes
473
       Module Procedure Fill_RefCodes_FAtom
474
       Module Procedure Fill_RefCodes_FmAtom
475
       Module Procedure Fill_RefCodes_Molcrys
476
       Module Procedure Fill_RefCodes_Molec
477
    End Interface
478

    
479
    Interface Get_ConCodes_Line
480
       Module Procedure Get_ConCodes_Line_FAtom
481
       Module Procedure Get_ConCodes_Line_FmAtom
482
       Module Procedure Get_ConCodes_Line_Molcrys
483
       Module Procedure Get_ConCodes_Line_Molec
484
    End Interface
485

    
486
    Interface Get_RefCodes_Line
487
       Module Procedure Get_RefCodes_Line_FAtom
488
       Module Procedure Get_RefCodes_Line_FmAtom
489
       Module Procedure Get_RefCodes_Line_Molcrys
490
       Module Procedure Get_RefCodes_Line_Molec
491
    End Interface
492

    
493
    Interface Init_RefCodes
494
       Module Procedure Init_RefCodes_FAtom
495
       Module Procedure Init_RefCodes_FmAtom
496
       Module Procedure Init_RefCodes_Molcrys
497
       Module Procedure Init_RefCodes_Molec
498
    End Interface
499

    
500
    Interface Read_RefCodes_File
501
       Module Procedure Read_RefCodes_File_FAtom
502
       Module Procedure Read_RefCodes_File_FmAtom
503
       Module Procedure Read_RefCodes_File_Molcrys
504
       Module Procedure Read_RefCodes_File_Molec
505
    End Interface
506

    
507
    Interface VState_to_AtomsPar
508
       Module Procedure VState_to_AtomsPar_FAtom
509
       Module Procedure VState_to_AtomsPar_FmAtom
510
       Module Procedure VState_to_AtomsPar_Molcrys
511
       Module Procedure VState_to_AtomsPar_Molec
512
    End Interface
513

    
514
    Interface Write_Info_RefCodes
515
       Module Procedure Write_Info_RefCodes_FAtom
516
       Module Procedure Write_Info_RefCodes_FmAtom
517
       Module Procedure Write_Info_RefCodes_Molcrys
518
       Module Procedure Write_Info_RefCodes_Molec
519
    End Interface
520

    
521
 Contains
522

    
523
    !!----
524
    !!---- Subroutine Allocate_RestParam(file_dat)
525
    !!----    Type(file_list_type),     intent( in)    :: file_dat
526
    !!----
527
    !!----    Allocate vectors Ang_Rest, Dist_Rest, Tor_Rest
528
    !!----
529
    !!---- Update: March - 2005
530
    !!
531
    Subroutine Allocate_RestParam(file_dat)
532
       !---- Arguments ----!
533
       Type(file_list_type),     intent( in)    :: file_dat
534

    
535
       !---- Local variables ----!
536
       character(len=132)              :: line
537
       character(len=15),dimension(40) :: car
538
       integer                         :: i,nc,nr
539

    
540
       if (allocated(Ang_Rest)) deallocate(Ang_Rest)
541
       if (allocated(Dis_Rest)) deallocate(Dis_Rest)
542
       if (allocated(Tor_Rest)) deallocate(Tor_Rest)
543

    
544
       NP_Rest_Ang=0
545
       NP_Rest_Dis=0
546
       NP_Rest_Tor=0
547

    
548
       !---- Dimension for AFIX ----!
549
       nr=0
550
       do i=1,file_dat%nlines
551
          line=adjustl(file_dat%line(i))
552
          if (u_case(line(1:4)) /= "AFIX") cycle
553
          call cutst(line)
554
          call getword(line,car,nc)
555
          nr=nr+nc/3
556
       end do
557
       if (nr >0) then
558
          allocate(Ang_Rest(nr))
559
          ang_rest%aobs =0.0
560
          ang_rest%acalc=0.0
561
          ang_rest%sigma=0.0
562
          ang_rest%p(1) = 0
563
          ang_rest%p(2) = 0
564
          ang_rest%STCode(1)=" "
565
          ang_rest%STCode(2)=" "
566
       end if
567

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

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

    
609
       return
610
    End Subroutine Allocate_RestParam
611

    
612
    !!---- Subroutine Allocate_VParam(N)
613
    !!----    integer, intent(in) :: N
614
    !!----
615
    !!----    Allocate vectors V_Vec, V_Bounds, V_Name, V_Bcon, V_Shift, V_list
616
    !!----    If N is equal zero it deallocates the vectors
617
    !!----
618
    !!---- Update: March - 2005
619
    !!
620
    Subroutine Allocate_VParam(N)
621
       !---- Arguments ----!
622
       integer, intent(in) :: N
623

    
624
       if (allocated(V_Vec))    deallocate(V_Vec)
625
       if (allocated(V_Name))   deallocate(V_Name)
626
       if (allocated(V_Bounds)) deallocate(V_Bounds)
627
       if (allocated(V_BCon))   deallocate(V_BCon)
628
       if (allocated(V_Shift))  deallocate(V_Shift)
629
       if (allocated(V_List))   deallocate(V_List)
630

    
631
       if (N > 0) then
632
          allocate(V_Vec(n))
633
          V_Vec=0.0
634
          allocate(V_Name(n))
635
          V_Name=" "
636
          allocate(V_Bounds(3,n))
637
          V_Bounds=0.0
638
          allocate(V_BCon(n))
639
          V_BCon=0
640
          allocate(V_Shift(n))
641
          V_Shift=0.0
642
          allocate(V_List(n))
643
          V_list=0
644
          np_max=n
645
       else
646
          np_max=0
647
       end if
648

    
649
       return
650
    End Subroutine Allocate_VParam
651

    
652
    !!--++
653
    !!--++ Subroutine Delete_RefCodes(N, FAtom/MolCrys/Molec)
654
    !!--++    integer,                      intent(in)     :: N
655
    !!--++    type(Atom_List_Type),         intent(in out) :: FAtom
656
    !!--++    or
657
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
658
    !!--++    or
659
    !!--++    type(molecule_type),          intent(in out) :: Molec
660
    !!--++
661
    !!--++    (Private)
662
    !!--++    Delete the number of Refinable Parameter (N) on the list
663
    !!--++
664
    !!--++ Update: March - 2005
665
    !!
666

    
667
    !!--++
668
    !!--++ Subroutine Delete_RefCodes_FAtom(N, FAtom)
669
    !!--++    integer,              intent(in)     :: N
670
    !!--++    type(Atom_List_Type), intent(in out) :: FAtom
671
    !!--++
672
    !!--++ Overloaded
673
    !!--++ Delete the number of Refinable Parameter (N) on the list
674
    !!--++
675
    !!--++ Update: March - 2005
676
    !!
677
    Subroutine Delete_RefCodes_FAtom(N, FAtom)
678
       !---- Arguments ----!
679
       integer,              intent(in)     :: N
680
       type(Atom_List_Type), intent(in out) :: FAtom
681

    
682
       !---- Local Variables ----!
683
       logical :: deleted
684
       integer :: i,j
685

    
686
       deleted=.false.
687

    
688
       !---- Eliminate N Parameter ----!
689
       do i=1,FAtom%natoms
690
          do j=1,3
691
             if (FAtom%atom(i)%lx(j) == N) then
692
                 FAtom%atom(i)%lx(j)=0
693
                 FAtom%atom(i)%mx(j)=0.0
694
                 deleted=.true.
695
             end if
696
          end do
697

    
698
          if (FAtom%atom(i)%lbiso == N) then
699
              FAtom%atom(i)%lbiso=0
700
              FAtom%atom(i)%mbiso=0.0
701
              deleted=.true.
702
          end if
703

    
704
          if (FAtom%atom(i)%locc == N) then
705
              FAtom%atom(i)%locc=0
706
              FAtom%atom(i)%mocc=0.0
707
              deleted=.true.
708
          end if
709

    
710
          do j=1,6
711
             if (FAtom%atom(i)%lu(j) == N) then
712
                 FAtom%atom(i)%lu(j)=0
713
                 FAtom%atom(i)%mu(j)=0.0
714
                 deleted=.true.
715
             end if
716
          end do
717
       end do
718

    
719
       !---- Updating Variables ----!
720
       do i=1,FAtom%natoms
721
          do j=1,3
722
             if (FAtom%atom(i)%lx(j) > N) then
723
                FAtom%atom(i)%lx(j)=FAtom%atom(i)%lx(j)-1
724
             end if
725
          end do
726

    
727
          if (FAtom%atom(i)%lbiso > N) then
728
             FAtom%atom(i)%lbiso=FAtom%atom(i)%lbiso-1
729
          end if
730

    
731
          if (FAtom%atom(i)%locc > N) then
732
             FAtom%atom(i)%locc=FAtom%atom(i)%locc-1
733
          end if
734

    
735
          do j=1,6
736
             if (FAtom%atom(i)%lu(j) > N) then
737
                FAtom%atom(i)%lu(j)=FAtom%atom(i)%lu(j)-1
738
             end if
739
          end do
740
       end do
741

    
742
       !---- Updating V_Vectors ----!
743
       if (deleted) then
744
          do i=N+1,Np_refi
745
             V_Vec(i-1)=V_Vec(i)
746
             V_Name(i-1)=V_Name(i)
747
             V_Bounds(:,i-1)=V_Bounds(:,i)
748
             V_BCon(i-1)=V_BCon(i)
749
             V_List(i-1)=V_List(i)
750
          end do
751
          V_Vec(np_refi)=0.0
752
          V_Name(np_refi)=" "
753
          V_Bounds(:,np_refi)=0.0
754
          V_BCon(np_refi)=0
755
          V_List(np_refi)=0
756

    
757
          np_refi=np_refi-1
758
       end if
759

    
760
       return
761
    End Subroutine Delete_RefCodes_FAtom
762

    
763
    !!--++
764
    !!--++ Subroutine Delete_RefCodes_FmAtom(N, FmAtom)
765
    !!--++    integer,              intent(in)     :: N
766
    !!--++    type(mAtom_List_Type),intent(in out) :: FmAtom
767
    !!--++
768
    !!--++ Delete the number of Refinable Parameters (N) on the list
769
    !!--++ magnetic clone of Delete_RefCodes_FAtom
770
    !!--++ Created: December - 2011
771
    !!
772
    Subroutine Delete_RefCodes_FmAtom(N, FmAtom)
773
       !---- Arguments ----!
774
       integer,              intent(in)     :: N
775
       type(mAtom_List_Type), intent(in out) :: FmAtom
776

    
777
       !---- Local Variables ----!
778
       logical :: deleted
779
       integer :: i,j,k,ik
780

    
781
       deleted=.false.
782

    
783
       !---- Eliminate N Parameter ----!
784
       do i=1,FmAtom%natoms
785
       ik=FmAtom%atom(i)%nvk
786
          do j=1,3
787
           do k=1,ik
788
             if (FmAtom%atom(i)%mSkR(j,k) == N) then
789
                 FmAtom%atom(i)%mSkR(j,k)=0.0
790
                 FmAtom%atom(i)%lskr(j,k)=0
791
                 deleted=.true.
792
             end if
793
           end do
794
          end do
795

    
796
          do j=1,3
797
           do k=1,ik
798
             if (FmAtom%atom(i)%mSkI(j,k) == N) then
799
                 FmAtom%atom(i)%mSkI(j,k)=0.0
800
                 FmAtom%atom(i)%lski(j,k)=0
801
                 deleted=.true.
802
             end if
803
           end do
804
          end do
805

    
806
          do k=1,ik
807
             if (FmAtom%atom(i)%mmphas(k) == N) then
808
                 FmAtom%atom(i)%mmphas(k)=0.0
809
                 FmAtom%atom(i)%lmphas(k)=0
810
                 deleted=.true.
811
             end if
812
          end do
813

    
814
          do j=1,12
815
           do k=1,ik
816
             if (FmAtom%atom(i)%mbas(j,k) == N) then
817
                 FmAtom%atom(i)%mbas(j,k)=0.0
818
                 FmAtom%atom(i)%lbas(j,k)=0
819
                 deleted=.true.
820
             end if
821
           end do
822
          end do
823
       end do
824

    
825
       !---- Updating Variables ----!
826
       do i=1,FmAtom%natoms
827
          do j=1,3
828
           do k=1,ik
829
             if (FmAtom%atom(i)%mSkR(j,k) > N) then
830
                 FmAtom%atom(i)%mSkR(j,k)=FmAtom%atom(i)%mSkR(j,k)-1
831
             end if
832
           end do
833
          end do
834
          
835
          do j=1,3
836
           do k=1,ik
837
             if (FmAtom%atom(i)%mSkI(j,k) > N) then
838
                 FmAtom%atom(i)%mSkI(j,k)=FmAtom%atom(i)%mSkI(j,k)-1
839
             end if
840
           end do
841
          end do
842

    
843
          do k=1,ik
844
             if (FmAtom%atom(i)%mmphas(k) > N) then
845
                 FmAtom%atom(i)%mmphas(k)=FmAtom%atom(i)%mmphas(k)-1
846
             end if
847
          end do
848

    
849
          do j=1,12
850
           do k=1,ik
851
             if (FmAtom%atom(i)%mbas(j,k) > N) then
852
                 FmAtom%atom(i)%mbas(j,k)=FmAtom%atom(i)%mbas(j,k)-1
853
             end if
854
           end do
855
          end do
856
       end do
857

    
858
       !---- Updating V_Vectors ----!
859
       if (deleted) then
860
          do i=N+1,Np_refi
861
             V_Vec(i-1)=V_Vec(i)
862
             V_Name(i-1)=V_Name(i)
863
             V_Bounds(:,i-1)=V_Bounds(:,i)
864
             V_BCon(i-1)=V_BCon(i)
865
             V_List(i-1)=V_List(i)
866
          end do
867
          V_Vec(np_refi)=0.0
868
          V_Name(np_refi)=" "
869
          V_Bounds(:,np_refi)=0.0
870
          V_BCon(np_refi)=0
871
          V_List(np_refi)=0
872

    
873
          np_refi=np_refi-1
874
       end if
875

    
876
       return
877
    End Subroutine Delete_RefCodes_FmAtom
878

    
879
    !!--++
880
    !!--++ Subroutine Delete_RefCodes_MolCrys(N,MolCrys)
881
    !!--++    integer,                      intent(in)     :: N
882
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
883
    !!--++
884
    !!--++ Overloaded
885
    !!--++ Delete the number of Refinable Parameter (N) on the list
886
    !!--++
887
    !!--++ Update: March - 2005
888
    !!
889
    Subroutine Delete_RefCodes_MolCrys(N,Molcrys)
890
       !---- Arguments ----!
891
       integer,                      intent(in)     :: N
892
       type(molecular_Crystal_type), intent(in out) :: MolCrys
893

    
894
       !---- Local Variables ----!
895
       logical :: deleted
896
       integer :: i,j,k
897

    
898
       deleted=.false.
899

    
900
       if (MolCrys%N_Free > 0 ) then
901
          do i=1,MolCrys%N_Free
902
             do j=1,3
903
                if (MolCrys%Atm(i)%lx(j) == N) then
904
                    MolCrys%Atm(i)%lx(j)=0
905
                    MolCrys%Atm(i)%mx(j)=0.0
906
                    deleted=.true.
907
                end if
908
             end do
909

    
910
             if (MolCrys%Atm(i)%lbiso == N) then
911
                 MolCrys%Atm(i)%lbiso=0
912
                 MolCrys%Atm(i)%mbiso=0.0
913
                 deleted=.true.
914
             end if
915

    
916
             if (MolCrys%Atm(i)%locc == N) then
917
                 MolCrys%Atm(i)%locc=0
918
                 MolCrys%Atm(i)%mocc=0.0
919
                 deleted=.true.
920
             end if
921

    
922
             do j=1,6
923
                if (MolCrys%Atm(i)%lu(j) == N) then
924
                    MolCrys%Atm(i)%lu(j)=0
925
                    MolCrys%Atm(i)%mu(j)=0.0
926
                    deleted=.true.
927
                end if
928
             end do
929
          end do
930

    
931
          do i=1,MolCrys%N_Free
932
             do j=1,3
933
                if (MolCrys%Atm(i)%lx(j) > N) then
934
                    MolCrys%Atm(i)%lx(j)=MolCrys%Atm(i)%lx(j)-1
935
                end if
936
             end do
937

    
938
             if (MolCrys%Atm(i)%lbiso > N) then
939
                MolCrys%Atm(i)%lbiso=MolCrys%Atm(i)%lbiso-1
940
             end if
941

    
942
             if (MolCrys%Atm(i)%locc > N) then
943
                MolCrys%Atm(i)%locc=MolCrys%Atm(i)%locc-1
944
             end if
945

    
946
             do j=1,6
947
                if (MolCrys%Atm(i)%lu(j) > N) then
948
                   MolCrys%Atm(i)%lu(j)=MolCrys%Atm(i)%lu(j)-1
949
                end if
950
             end do
951
          end do
952
       end if
953

    
954
       if (MolCrys%N_Mol > 0 ) then
955

    
956
          do k=1,MolCrys%N_Mol
957
             do j=1,3
958
                if (Molcrys%Mol(k)%lxcentre(j) == N) then
959
                   Molcrys%Mol(k)%lxcentre(j)=0
960
                   Molcrys%Mol(k)%mxcentre(j)=0.0
961
                   deleted=.true.
962
                end if
963
             end do
964

    
965
             do j=1,3
966
                if (Molcrys%Mol(k)%lOrient(j) == N) then
967
                   Molcrys%Mol(k)%lOrient(j)=0
968
                   Molcrys%Mol(k)%mOrient(j)=0.0
969
                   deleted=.true.
970
                end if
971
             end do
972

    
973
             do j=1,6
974
                if (Molcrys%Mol(k)%lT_TLS(j) == N) then
975
                   Molcrys%Mol(k)%lT_TLS(j)=0
976
                   Molcrys%Mol(k)%mT_TLS(j)=0.0
977
                   deleted=.true.
978
                end if
979
             end do
980

    
981
             do j=1,6
982
                if (Molcrys%Mol(k)%lL_TLS(j) == N) then
983
                   Molcrys%Mol(k)%lL_TLS(j)=0
984
                   Molcrys%Mol(k)%mL_TLS(j)=0.0
985
                   deleted=.true.
986
                end if
987
             end do
988

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

    
999
             !---- Updating ----!
1000
             do j=1,3
1001
                if (Molcrys%Mol(k)%lxcentre(j) > N) then
1002
                   Molcrys%Mol(k)%lxcentre(j)=Molcrys%Mol(k)%lxcentre(j)-1
1003
                end if
1004
             end do
1005

    
1006
             do j=1,3
1007
                if (Molcrys%Mol(k)%lOrient(j) > N) then
1008
                   Molcrys%Mol(k)%lOrient(j)=Molcrys%Mol(k)%lOrient(j)-1
1009
                end if
1010
             end do
1011

    
1012
             do j=1,6
1013
                if (Molcrys%Mol(k)%lT_TLS(j) > N) then
1014
                   Molcrys%Mol(k)%lT_TLS(j)=Molcrys%Mol(k)%lT_TLS(j)-1
1015
                end if
1016
             end do
1017

    
1018
             do j=1,6
1019
                if (Molcrys%Mol(k)%lL_TLS(j) > N) then
1020
                   Molcrys%Mol(k)%lL_TLS(j)=Molcrys%Mol(k)%lL_TLS(j)-1
1021
                end if
1022
             end do
1023

    
1024
             do i=1,3
1025
                do j=1,3
1026
                   if (Molcrys%Mol(k)%lS_TLS(i,j) > N) then
1027
                      Molcrys%Mol(k)%lS_TLS(i,j)=Molcrys%Mol(k)%lS_TLS(i,j)-1
1028
                   end if
1029
                end do
1030
             end do
1031

    
1032
             if (Molcrys%Mol(k)%natoms <=0) cycle
1033

    
1034
             do i=1,Molcrys%Mol(k)%natoms
1035
                do j=1,3
1036
                   if (MolCrys%Mol(k)%lI_coor(j,i) == N) then
1037
                      MolCrys%Mol(k)%lI_coor(j,i)=0
1038
                      MolCrys%Mol(k)%mI_coor(j,i)=0.0
1039
                      deleted=.true.
1040
                   end if
1041
                end do
1042

    
1043
                if (MolCrys%Mol(k)%lbiso(i) == N) then
1044
                   MolCrys%Mol(k)%lbiso(i)=0
1045
                   MolCrys%Mol(k)%mbiso(i)=0.0
1046
                   deleted=.true.
1047
                end if
1048

    
1049
                if (MolCrys%Mol(k)%locc(i) == N) then
1050
                   MolCrys%Mol(k)%locc(i)=0
1051
                   MolCrys%Mol(k)%mocc(i)=0.0
1052
                   deleted=.true.
1053
                end if
1054
             end do
1055

    
1056
             do i=1,Molcrys%Mol(k)%natoms
1057
                do j=1,3
1058
                   if (MolCrys%Mol(k)%lI_coor(j,i) > N) then
1059
                      MolCrys%Mol(k)%lI_coor(j,i)=MolCrys%Mol(k)%lI_coor(j,i)-1
1060
                   end if
1061
                end do
1062

    
1063
                if (MolCrys%Mol(k)%lbiso(i) > N) then
1064
                   MolCrys%Mol(k)%lbiso(i)=MolCrys%Mol(k)%lbiso(i)-1
1065
                end if
1066

    
1067
                if (MolCrys%Mol(k)%locc(i) > N) then
1068
                   MolCrys%Mol(k)%locc(i)=MolCrys%Mol(k)%locc(i)-1
1069
                end if
1070
             end do
1071

    
1072
          end do
1073
       end if
1074

    
1075
       !---- Updating V_Vectors ----!
1076
       if (deleted) then
1077
          do i=N+1,Np_refi
1078
             V_Vec(i-1)=V_Vec(i)
1079
             V_Name(i-1)=V_Name(i)
1080
             V_Bounds(:,i-1)=V_Bounds(:,i)
1081
             V_BCon(i-1)=V_BCon(i)
1082
             V_List(i-1)=V_List(i)
1083
          end do
1084
          V_Vec(np_refi)=0.0
1085
          V_Name(np_refi)=" "
1086
          V_Bounds(:,np_refi)=0.0
1087
          V_BCon(np_refi)=0
1088
          V_List(np_refi)=0
1089

    
1090
          np_refi=np_refi-1
1091
       end if
1092

    
1093
       return
1094
    End Subroutine Delete_RefCodes_MolCrys
1095

    
1096
    !!--++
1097
    !!--++ Subroutine Delete_RefCodes_Molec(N,Molec)
1098
    !!--++    integer,             intent(in)     :: N
1099
    !!--++    type(molecule_type), intent(in out) :: Molec
1100
    !!--++
1101
    !!--++ Overloaded
1102
    !!--++ Delete the number of Refinable Parameter (N) on the list
1103
    !!--++
1104
    !!--++ Update: March - 2005
1105
    !!
1106
    Subroutine Delete_RefCodes_Molec(N,Molec)
1107
       !---- Arguments ----!
1108
       integer,             intent(in)     :: N
1109
       type(molecule_type), intent(in out) :: Molec
1110

    
1111
       !---- Local Variables ----!
1112
       logical :: deleted
1113
       integer :: i,j
1114

    
1115
       deleted=.false.
1116

    
1117
       do j=1,3
1118
          if (Molec%lxcentre(j) == N) then
1119
             Molec%lxcentre(j)=0
1120
             Molec%mxcentre(j)=0.0
1121
             deleted=.true.
1122
          end if
1123
       end do
1124

    
1125
       do j=1,3
1126
          if (Molec%lOrient(j) == N) then
1127
             Molec%lOrient(j)=0
1128
             Molec%mOrient(j)=0.0
1129
             deleted=.true.
1130
          end if
1131
       end do
1132

    
1133
       do j=1,6
1134
          if (Molec%lT_TLS(j) == N) then
1135
             Molec%lT_TLS(j)=0
1136
             Molec%mT_TLS(j)=0.0
1137
             deleted=.true.
1138
          end if
1139
       end do
1140

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

    
1149
       do i=1,3
1150
          do j=1,3
1151
             if (Molec%lS_TLS(i,j) == N) then
1152
                Molec%lS_TLS(i,j)=0
1153
                Molec%mS_TLS(i,j)=0.0
1154
                deleted=.true.
1155
             end if
1156
          end do
1157
       end do
1158

    
1159
       !---- Updating ----!
1160
       do j=1,3
1161
          if (Molec%lxcentre(j) > N) then
1162
             Molec%lxcentre(j)=Molec%lxcentre(j)-1
1163
          end if
1164
       end do
1165

    
1166
       do j=1,3
1167
          if (Molec%lOrient(j) > N) then
1168
             Molec%lOrient(j)=Molec%lOrient(j)-1
1169
          end if
1170
       end do
1171

    
1172
       do j=1,6
1173
          if (Molec%lT_TLS(j) > N) then
1174
             Molec%lT_TLS(j)=Molec%lT_TLS(j)-1
1175
          end if
1176
       end do
1177

    
1178
       do j=1,6
1179
          if (Molec%lL_TLS(j) > N) then
1180
             Molec%lL_TLS(j)=Molec%lL_TLS(j)-1
1181
          end if
1182
       end do
1183

    
1184
       do i=1,3
1185
          do j=1,3
1186
             if (Molec%lS_TLS(i,j) > N) then
1187
                Molec%lS_TLS(i,j)=Molec%lS_TLS(i,j)-1
1188
             end if
1189
          end do
1190
       end do
1191

    
1192
       if (molec%natoms <=0) return
1193

    
1194
       do i=1,Molec%Natoms
1195
          do j=1,3
1196
             if (Molec%lI_coor(j,i) == N) then
1197
                Molec%lI_coor(j,i)=0
1198
                Molec%mI_coor(j,i)=0.0
1199
                deleted=.true.
1200
             end if
1201
          end do
1202

    
1203
          if (Molec%lbiso(i) == N) then
1204
             Molec%lbiso(i)=0
1205
             Molec%mbiso(i)=0.0
1206
             deleted=.true.
1207
          end if
1208

    
1209
          if (Molec%locc(i) == N) then
1210
             Molec%locc(i)=0
1211
             Molec%mocc(i)=0.0
1212
             deleted=.true.
1213
          end if
1214
       end do
1215

    
1216
       !---- Updating ----!
1217
       do i=1,Molec%Natoms
1218
          do j=1,3
1219
             if (Molec%lI_coor(j,i) > N) then
1220
                Molec%lI_coor(j,i)=Molec%lI_coor(j,i)-1
1221
             end if
1222
          end do
1223

    
1224
          if (Molec%lbiso(i) > N) then
1225
             Molec%lbiso(i)=Molec%lbiso(i)-1
1226
          end if
1227

    
1228
          if (Molec%locc(i) > N) then
1229
             Molec%locc(i)=Molec%locc(i)-1
1230
          end if
1231
       end do
1232

    
1233
       !---- Updating V_Vectors ----!
1234
       if (deleted) then
1235
          do i=N+1,Np_refi
1236
             V_Vec(i-1)=V_Vec(i)
1237
             V_Name(i-1)=V_Name(i)
1238
             V_Bounds(:,i-1)=V_Bounds(:,i)
1239
             V_BCon(i-1)=V_BCon(i)
1240
             V_List(i-1)=V_List(i)
1241
          end do
1242
          V_Vec(np_refi)=0.0
1243
          V_Name(np_refi)=" "
1244
          V_Bounds(:,np_refi)=0.0
1245
          V_BCon(np_refi)=0
1246
          V_List(np_refi)=0
1247

    
1248
          np_refi=np_refi-1
1249
       end if
1250

    
1251
       return
1252
    End Subroutine Delete_RefCodes_Molec
1253

    
1254
    !!--++
1255
    !!--++ Subroutine Fill_RefCodes(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom/MolCrys/Molec,Spg)
1256
    !!--++    integer,                       intent(in)     :: Key
1257
    !!--++    character(len=*),              intent(in)     :: Dire
1258
    !!--++    integer,                       intent(in)     :: Na
1259
    !!--++    integer,                       intent(in)     :: Nb
1260
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1261
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1262
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1263
    !!--++    integer,                       intent(in)     :: Ic
1264
    !!--++    type(Atom_List_Type),          intent(in out) :: FAtom
1265
    !!--++    or
1266
    !!--++    type(molecular_Crystal_type),  intent(in out) :: MolCrys
1267
    !!--++    or
1268
    !!--++    type(molecule_type),           intent(in out) :: Molec
1269
    !!--++    type(space_group_type),        intent(in)     :: Spg
1270
    !!--++
1271
    !!--++    (Private)
1272
    !!--++    Write on Vectors the Information for Free Atoms
1273
    !!--++
1274
    !!--++ Update: March - 2005
1275
    !!
1276

    
1277
    !!--++
1278
    !!--++ Subroutine Fill_RefCodes_FAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom,Spg)
1279
    !!--++    integer,                       intent(in)     :: Key
1280
    !!--++    character(len=*),              intent(in)     :: Dire
1281
    !!--++    integer,                       intent(in)     :: Na
1282
    !!--++    integer,                       intent(in)     :: Nb
1283
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1284
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1285
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1286
    !!--++    integer,                       intent(in)     :: Ic
1287
    !!--++    type(Atom_List_Type),          intent(in out) :: FAtom
1288
    !!--++    type(space_group_type),        intent(in)     :: Spg
1289
    !!--++
1290
    !!--++ Overloaded
1291
    !!--++ Write on Vectors the Information for Free Atoms
1292
    !!--++
1293
    !!--++ Update: March - 2005
1294
    !!
1295
    Subroutine Fill_RefCodes_FAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FAtom,Spg)
1296
       !---- Arguments ----!
1297
       integer,                       intent(in)     :: Key
1298
       character(len=*),              intent(in)     :: Dire
1299
       integer,                       intent(in)     :: Na
1300
       integer,                       intent(in)     :: Nb
1301
       real(kind=cp),                 intent(in)     :: Xl
1302
       real(kind=cp),                 intent(in)     :: Xu
1303
       real(kind=cp),                 intent(in)     :: Xs
1304
       integer,                       intent(in)     :: Ic
1305
       type(Atom_List_Type),          intent(in out) :: FAtom
1306
       type(space_group_type),        intent(in)     :: Spg
1307

    
1308
       !---- Local variables ----!
1309
       integer           :: j,nc
1310

    
1311
       call init_err_refcodes()
1312
       if (Na <= 0) then
1313
          err_refcodes=.true.
1314
          ERR_RefCodes_Mess="Number of atom no defined"
1315
          return
1316
       end if
1317

    
1318
       select case (dire)
1319
          !---- FIX Directive ----!
1320
          case ("fix")
1321

    
1322
             select case (key)
1323
                case (0)
1324
                   !---- nb must be different zero ----!
1325
                   select case (nb)
1326
                      case (0)
1327
                         err_refcodes=.true.
1328
                         ERR_RefCodes_Mess="Option not defined"
1329
                         return
1330

    
1331
                      case ( 1:3)
1332
                         !---- X_, Y_, Z_ ----!
1333
                         if (FAtom%atom(na)%lx(nb) /=0) then
1334
                            nc=FAtom%atom(na)%lx(nb)
1335
                            call Delete_RefCodes(nc,FAtom)
1336
                         end if
1337

    
1338
                      case ( 4)
1339
                         !---- Biso_ ----!
1340
                         if (FAtom%atom(na)%lbiso /=0) then
1341
                            nc=FAtom%atom(na)%lbiso
1342
                            call Delete_RefCodes(nc,fatom)
1343
                         end if
1344

    
1345
                      case ( 5)
1346
                         !---- Occ_ ----!
1347
                         if (FAtom%atom(na)%locc /=0) then
1348
                            nc=FAtom%atom(na)%locc
1349
                            call Delete_RefCodes(nc,fatom)
1350
                         end if
1351

    
1352
                      case ( 6:11)
1353
                         !---- B11_,...,B23_ ----!
1354
                         if (FAtom%atom(na)%lu(nb-5) /=0) then
1355
                            nc=FAtom%atom(na)%lu(nb-5)
1356
                            call Delete_RefCodes(nc,fatom)
1357
                         end if
1358

    
1359
                      case (12)
1360
                         !---- Banis_ ----!
1361
                         do j=1,6
1362
                            if (FAtom%atom(na)%lu(j) /=0) then
1363
                               nc=FAtom%atom(na)%lu(j)
1364
                               call Delete_RefCodes(nc,fatom)
1365
                            end if
1366
                         end do
1367

    
1368
                      case (13:)
1369
                         err_refcodes=.true.
1370
                         ERR_RefCodes_Mess="Option not defined for this type of variable "
1371
                         return
1372
                   end select ! nb
1373

    
1374
                case (1)
1375
                   !---- XYZ ----!
1376
                   do j=1,3
1377
                      if (FAtom%atom(na)%lx(j) /=0) then
1378
                         nc=FAtom%atom(na)%lx(j)
1379
                         call Delete_RefCodes(nc,fatom)
1380
                      end if
1381
                   end do
1382

    
1383
                case (2)
1384
                   !---- OCC ----!
1385
                   if (FAtom%atom(na)%locc /=0) then
1386
                      nc=FAtom%atom(na)%locc
1387
                      call Delete_RefCodes(nc,fatom)
1388
                   end if
1389

    
1390
                case (3)
1391
                   !---- BIS ----!
1392
                   if (FAtom%atom(na)%lbiso /=0) then
1393
                      nc=FAtom%atom(na)%lbiso
1394
                      call Delete_RefCodes(nc,fatom)
1395
                   end if
1396

    
1397
                case (4)
1398
                  !---- BAN ----!
1399
                  do j=1,6
1400
                     if (FAtom%atom(na)%lu(j) /=0) then
1401
                        nc=FAtom%atom(na)%lu(j)
1402
                        call Delete_RefCodes(nc,fatom)
1403
                     end if
1404
                  end do
1405

    
1406
                case (5)
1407
                   !---- ALL ----!
1408
                   do j=1,3
1409
                      if (FAtom%atom(na)%lx(j) /=0) then
1410
                         nc=FAtom%atom(na)%lx(j)
1411
                         call Delete_RefCodes(nc,fatom)
1412
                      end if
1413
                   end do
1414
                   if (FAtom%atom(na)%locc /=0) then
1415
                      nc=FAtom%atom(na)%locc
1416
                      call Delete_RefCodes(nc,fatom)
1417
                   end if
1418
                   if (FAtom%atom(na)%lbiso /=0) then
1419
                      nc=FAtom%atom(na)%lbiso
1420
                      call Delete_RefCodes(nc,fatom)
1421
                   end if
1422
                   do j=1,6
1423
                      if (FAtom%atom(na)%lu(j) /=0) then
1424
                         nc=FAtom%atom(na)%lu(j)
1425
                         call Delete_RefCodes(nc,fatom)
1426
                      end if
1427
                   end do
1428

    
1429
                case (6:)
1430
                   err_refcodes=.true.
1431
                   ERR_RefCodes_Mess="Incompatible information for this type of variable "
1432
                   return
1433
             end select
1434

    
1435
          !---- VARY Directive ----!
1436
          case ("var")
1437

    
1438
             select case (key)
1439
                case (0)
1440

    
1441
                   !---- nb must be different zero ----!
1442
                   select case (nb)
1443
                      case (0)
1444
                         err_refcodes=.true.
1445
                         ERR_RefCodes_Mess="Option not defined"
1446
                         return
1447

    
1448
                      case ( 1:3)
1449
                         !---- X_, Y_, Z_ ----!
1450
                         if (FAtom%atom(na)%lx(nb) ==0) then
1451
                            FAtom%atom(na)%mx(nb)=1.0
1452
                            call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1453
                                                 FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1454
                            if (FAtom%atom(na)%lx(nb) == np_refi) then
1455
                               V_Vec(np_refi)=FAtom%atom(na)%x(nb)
1456
                               V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1457
                               V_Bounds(1,np_refi)=xl
1458
                               V_Bounds(2,np_refi)=xu
1459
                               V_Bounds(3,np_refi)=xs
1460
                               V_BCon(np_refi)=ic
1461
                               V_list(np_refi)=na
1462
                            else
1463
                               np_refi=np_refi-1
1464
                            end if
1465
                         end if
1466

    
1467
                      case ( 4)
1468
                         !---- Biso_ ----!
1469
                         if (FAtom%atom(na)%lbiso ==0) then
1470
                            np_refi=np_refi+1
1471
                            V_Vec(np_refi)=FAtom%atom(na)%biso
1472
                            V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1473
                            FAtom%atom(na)%mbiso=1.0
1474
                            FAtom%atom(na)%lbiso=np_refi
1475
                            V_Bounds(1,np_refi)=xl
1476
                            V_Bounds(2,np_refi)=xu
1477
                            V_Bounds(3,np_refi)=xs
1478
                            V_BCon(np_refi)=ic
1479
                            V_list(np_refi)=na
1480
                         end if
1481

    
1482
                      case ( 5)
1483
                         !---- Occ_ ----!
1484
                         if (FAtom%atom(na)%locc ==0) then
1485
                            np_refi=np_refi+1
1486
                            V_Vec(np_refi)=FAtom%atom(na)%occ
1487
                            V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1488
                            FAtom%atom(na)%mocc=1.0
1489
                            FAtom%atom(na)%locc=np_refi
1490
                            V_Bounds(1,np_refi)=xl
1491
                            V_Bounds(2,np_refi)=xu
1492
                            V_Bounds(3,np_refi)=xs
1493
                            V_BCon(np_refi)=ic
1494
                            V_list(np_refi)=na
1495
                         end if
1496

    
1497
                      case ( 6:11)
1498
                         !---- B11_,...,B23_ ----!
1499
                         if (FAtom%atom(na)%lu(nb-5) ==0) then
1500
                            FAtom%atom(na)%mu(nb-5)=1.0
1501
                            call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1502
                                                 np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1503
                            if (FAtom%atom(na)%lu(nb-5) == np_refi) then
1504
                               V_Vec(np_refi)=FAtom%atom(na)%u(nb-5)
1505
                               V_Name(np_refi)=trim(code_nam(nb))//trim(FAtom%atom(na)%lab)
1506
                               V_Bounds(1,np_refi)=xl
1507
                               V_Bounds(2,np_refi)=xu
1508
                               V_Bounds(3,np_refi)=xs
1509
                               V_BCon(np_refi)=ic
1510
                               V_list(np_refi)=na
1511
                            else
1512
                               np_refi=np_refi-1
1513
                            end if
1514
                         end if
1515

    
1516
                      case (12)
1517
                         !---- Banis_ ----!
1518
                         do j=1,6
1519
                            if (FAtom%atom(na)%lu(j) ==0) then
1520
                               FAtom%atom(na)%mu(j)=1.0
1521
                               call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1522
                                                    np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1523
                               if (FAtom%atom(na)%lu(j) == np_refi) then
1524
                                  V_Vec(np_refi)=FAtom%atom(na)%u(j)
1525
                                  V_Name(np_refi)=trim(code_nam(5+j))//trim(FAtom%atom(na)%lab)
1526
                                  V_Bounds(1,np_refi)=xl
1527
                                  V_Bounds(2,np_refi)=xu
1528
                                  V_Bounds(3,np_refi)=xs
1529
                                  V_BCon(np_refi)=ic
1530
                                  V_list(np_refi)=na
1531
                               else
1532
                                  np_refi=np_refi-1
1533
                               end if
1534
                            end if
1535
                         end do
1536

    
1537
                      case (13:)
1538
                         err_refcodes=.true.
1539
                         ERR_RefCodes_Mess="Option Not defined by this type of variables"
1540
                         return
1541

    
1542
                   end select ! nb
1543

    
1544
                case (1)
1545
                   !---- XYZ ----!
1546
                   do j=1,3
1547
                      if (FAtom%atom(na)%lx(j) ==0) then
1548
                         FAtom%atom(na)%mx(j)=1.0
1549
                         call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1550
                                              FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1551
                         if (FAtom%atom(na)%lx(j) == np_refi) then
1552
                            V_Vec(np_refi)=FAtom%atom(na)%x(j)
1553
                            V_Name(np_refi)=trim(code_nam(j))//trim(FAtom%atom(na)%lab)
1554
                            V_Bounds(1,np_refi)=xl
1555
                            V_Bounds(2,np_refi)=xu
1556
                            V_Bounds(3,np_refi)=xs
1557
                            V_BCon(np_refi)=ic
1558
                            V_list(np_refi)=na
1559
                         else
1560
                            np_refi=np_refi-1
1561
                         end if
1562
                      end if
1563
                   end do
1564

    
1565
                case (2)
1566
                   !---- OCC ----!
1567
                   if (FAtom%atom(na)%locc ==0) then
1568
                      np_refi=np_refi+1
1569
                      V_Vec(np_refi)=FAtom%atom(na)%occ
1570
                      V_Name(np_refi)=trim(code_nam(5))//trim(FAtom%atom(na)%lab)
1571
                      FAtom%atom(na)%mocc=1.0
1572
                      FAtom%atom(na)%locc=np_refi
1573
                      V_Bounds(1,np_refi)=xl
1574
                      V_Bounds(2,np_refi)=xu
1575
                      V_Bounds(3,np_refi)=xs
1576
                      V_BCon(np_refi)=ic
1577
                      V_list(np_refi)=na
1578
                   end if
1579

    
1580
                case (3)
1581
                   !---- BIS ----!
1582
                   if (FAtom%atom(na)%lbiso ==0) then
1583
                      np_refi=np_refi+1
1584
                      V_Vec(np_refi)=FAtom%atom(na)%biso
1585
                      V_Name(np_refi)=trim(code_nam(4))//trim(FAtom%atom(na)%lab)
1586
                      FAtom%atom(na)%mbiso=1.0
1587
                      FAtom%atom(na)%lbiso=np_refi
1588
                      V_Bounds(1,np_refi)=xl
1589
                      V_Bounds(2,np_refi)=xu
1590
                      V_Bounds(3,np_refi)=xs
1591
                      V_BCon(np_refi)=ic
1592
                      V_list(np_refi)=na
1593
                   end if
1594

    
1595
                case (4)
1596
                   !---- BAN ----!
1597
                   do j=1,6
1598
                      if (FAtom%atom(na)%lu(j) ==0) then
1599
                         FAtom%atom(na)%mu(j)=1.0
1600
                         call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1601
                                              np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1602
                         if (FAtom%atom(na)%lu(j) == np_refi) then
1603
                            V_Vec(np_refi)=FAtom%atom(na)%u(j)
1604
                            V_Name(np_refi)=trim(code_nam(5+j))//trim(FAtom%atom(na)%lab)
1605
                            V_Bounds(1,np_refi)=xl
1606
                            V_Bounds(2,np_refi)=xu
1607
                            V_Bounds(3,np_refi)=xs
1608
                            V_BCon(np_refi)=ic
1609
                            V_list(np_refi)=na
1610
                         else
1611
                            np_refi=np_refi-1
1612
                         end if
1613
                      end if
1614
                   end do
1615

    
1616
                case (5)
1617
                   !---- ALL ----!
1618
                   do j=1,3
1619
                      if (FAtom%atom(na)%lx(j) ==0) then
1620
                         FAtom%atom(na)%mx(j)=1.0
1621
                         call get_atompos_ctr(FAtom%atom(na)%x, Spg, np_refi,   &
1622
                                              FAtom%atom(na)%lx, FAtom%atom(na)%mx)
1623
                         if (FAtom%atom(na)%lx(j) == np_refi) then
1624
                            V_Vec(np_refi)=FAtom%atom(na)%x(j)
1625
                            V_Name(np_refi)=trim(code_nam(j))//trim(FAtom%atom(na)%lab)
1626
                            V_Bounds(1,np_refi)=xl
1627
                            V_Bounds(2,np_refi)=xu
1628
                            V_Bounds(3,np_refi)=xs
1629
                            V_BCon(np_refi)=ic
1630
                            V_list(np_refi)=na
1631
                         else
1632
                            np_refi=np_refi-1
1633
                         end if
1634
                      end if
1635
                   end do
1636
                   if (FAtom%atom(na)%locc ==0) then
1637
                      np_refi=np_refi+1
1638
                      V_Vec(np_refi)=FAtom%atom(na)%occ
1639
                      V_Name(np_refi)=trim(mcode_nam(5))//trim(FAtom%atom(na)%lab)
1640
                      FAtom%atom(na)%mocc=1.0
1641
                      FAtom%atom(na)%locc=np_refi
1642
                      V_Bounds(1,np_refi)=xl
1643
                      V_Bounds(2,np_refi)=xu
1644
                      V_Bounds(3,np_refi)=xs
1645
                      V_BCon(np_refi)=ic
1646
                      V_list(np_refi)=na
1647
                   end if
1648
                   if (FAtom%atom(na)%lbiso ==0) then
1649
                      np_refi=np_refi+1
1650
                      V_Vec(np_refi)=FAtom%atom(na)%biso
1651
                      V_Name(np_refi)=trim(mcode_nam(4))//trim(FAtom%atom(na)%lab)
1652
                      FAtom%atom(na)%mbiso=1.0
1653
                      FAtom%atom(na)%lbiso=np_refi
1654
                      V_Bounds(1,np_refi)=xl
1655
                      V_Bounds(2,np_refi)=xu
1656
                      V_Bounds(3,np_refi)=xs
1657
                      V_BCon(np_refi)=ic
1658
                      V_list(np_refi)=na
1659
                   end if
1660
                   do j=1,6
1661
                      if (FAtom%atom(na)%lu(j) ==0) then
1662
                         FAtom%atom(na)%mu(j)=1.0
1663
                         call get_atombet_ctr(FAtom%atom(na)%x,FAtom%atom(na)%u,Spg, &
1664
                                              np_refi,FAtom%atom(na)%lu,FAtom%atom(na)%mu)
1665
                         if (FAtom%atom(na)%lu(j) == np_refi) then
1666
                            V_Vec(np_refi)=FAtom%atom(na)%u(j)
1667
                            V_Name(np_refi)=trim(mcode_nam(5+j))//trim(FAtom%atom(na)%lab)
1668
                            V_Bounds(1,np_refi)=xl
1669
                            V_Bounds(2,np_refi)=xu
1670
                            V_Bounds(3,np_refi)=xs
1671
                            V_BCon(np_refi)=ic
1672
                            V_list(np_refi)=na
1673
                         else
1674
                            np_refi=np_refi-1
1675
                         end if
1676
                      end if
1677
                   end do
1678

    
1679
                case(6:)
1680
                   err_refcodes=.true.
1681
                   ERR_RefCodes_Mess="Option Not defined by this type of variables"
1682
                   return
1683
             end select
1684
       end select
1685

    
1686
       return
1687
    End Subroutine Fill_RefCodes_FAtom
1688

    
1689
    !!--++
1690
    !!--++ Subroutine Fill_RefCodes_FmAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FmAtom,Spg)
1691
    !!--++    integer,                       intent(in)     :: Key
1692
    !!--++    character(len=*),              intent(in)     :: Dire
1693
    !!--++    integer,                       intent(in)     :: Na
1694
    !!--++    integer,                       intent(in)     :: Nb
1695
    !!--++    real(kind=cp),                 intent(in)     :: Xl
1696
    !!--++    real(kind=cp),                 intent(in)     :: Xu
1697
    !!--++    real(kind=cp),                 intent(in)     :: Xs
1698
    !!--++    integer,                       intent(in)     :: Ic
1699
    !!--++    type(mAtom_List_Type),         intent(in out) :: FmAtom
1700
    !!--++    type(space_group_type),        intent(in)     :: Spg
1701
    !!--++
1702
    !!--++ Write on Vectors the Information for Free Magnetic Atoms
1703
    !!--++ magnetic clone of Fill_RefCodes_FAtom
1704
    !!--++ Created: December - 2011
1705
    !!--++
1706
    Subroutine Fill_RefCodes_FmAtom(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,FmAtom,Spg)
1707
       !---- Arguments ----!
1708
       integer,                       intent(in)     :: Key
1709
       character(len=*),              intent(in)     :: Dire
1710
       integer,                       intent(in)     :: Na
1711
       integer,                       intent(in)     :: Nb
1712
       real(kind=cp),                 intent(in)     :: Xl
1713
       real(kind=cp),                 intent(in)     :: Xu
1714
       real(kind=cp),                 intent(in)     :: Xs
1715
       integer,                       intent(in)     :: Ic
1716
       type(mAtom_List_Type),         intent(in out) :: FmAtom
1717
       type(space_group_type),        intent(in)     :: Spg
1718

    
1719
       !---- Local variables ----!
1720
       integer           :: j,nc,ik
1721

    
1722
       call init_err_refcodes()
1723
       if (Na <= 0) then
1724
          err_refcodes=.true.
1725
          ERR_RefCodes_Mess="Number of atom no defined"
1726
          return
1727
       end if
1728

    
1729
       !---- Get im, number of the magnetic matrices/irrep set
1730
       ik=FmAtom%atom(na)%nvk
1731

    
1732
       select case (dire)
1733
          !---- FIX Directive ----!
1734
          case ("fix")
1735

    
1736
             select case (key)
1737
                case (0)
1738
                   !---- nb must be different zero ----!
1739
                   select case (nb)
1740
                      case (0)
1741
                         err_refcodes=.true.
1742
                         ERR_RefCodes_Mess="Option not defined"
1743
                         return
1744

    
1745
                      case ( 1:3)
1746
                         !---- Rx_, Ry_, Rz_ ----!
1747
                         if (FmAtom%atom(na)%lSkR(nb,ik) /=0) then
1748
                            nc=FmAtom%atom(na)%lSkR(nb,ik)
1749
                            call Delete_RefCodes(nc,FmAtom)
1750
                         end if
1751

    
1752
                      case ( 4:6)
1753
                         !---- Ix_, Iy_, Iz_ ----!
1754
                         if (FmAtom%atom(na)%lSkI(nb-3,ik) /=0) then
1755
                            nc=FmAtom%atom(na)%lSkI(nb-3,ik)
1756
                            call Delete_RefCodes(nc,FmAtom)
1757
                         end if
1758

    
1759
                      case ( 7:9)
1760
                         !---- Rm_, Rphi_, Rth_ ----!
1761
                         if (FmAtom%atom(na)%lSkR(nb-6,ik) /=0) then
1762
                            nc=FmAtom%atom(na)%lSkR(nb-6,ik)
1763
                            call Delete_RefCodes(nc,FmAtom)
1764
                         end if
1765

    
1766
                      case (10:12)
1767
                         !---- Im_, Iphi_, Ith_ ----!
1768
                         if (FmAtom%atom(na)%lSkI(nb-9,ik) /=0) then
1769
                            nc=FmAtom%atom(na)%lSkI(nb-9,ik)
1770
                            call Delete_RefCodes(nc,FmAtom)
1771
                         end if
1772

    
1773
                      case (13)
1774
                         !---- MagPh_ ----!
1775
                            if (FmAtom%atom(na)%lmphas(ik) /=0) then
1776
                               nc=FmAtom%atom(na)%lmphas(ik)
1777
                               call Delete_RefCodes(nc,FmAtom)
1778
                            end if
1779

    
1780
                      case (14:22)
1781
                         !---- C1_,..C12_ ----!
1782
                            if (FmAtom%atom(na)%lbas(nb-13,ik) /=0) then
1783
                               nc=FmAtom%atom(na)%lbas(nb-13,ik)
1784
                               call Delete_RefCodes(nc,FmAtom)
1785
                            end if
1786

    
1787
                      case (23)
1788
                         err_refcodes=.true.
1789
                         ERR_RefCodes_Mess="Option not defined for this type of variable "
1790
                         return
1791
                   end select ! nb
1792

    
1793
                case (1)
1794
                   !---- Rxyz ----!
1795
                   do j=1,3
1796
                      if (FmAtom%atom(na)%lSkR(j,ik) /=0) then
1797
                         nc=FmAtom%atom(na)%lSkR(j,ik)
1798
                         call Delete_RefCodes(nc,FmAtom)
1799
                      end if
1800
                   end do
1801

    
1802
                case (2)
1803
                   !---- Ixyz ----!
1804
                   do j=1,3
1805
                      if (FmAtom%atom(na)%lSkI(j,ik) /=0) then
1806
                         nc=FmAtom%atom(na)%lSkI(j,ik)
1807
                         call Delete_RefCodes(nc,FmAtom)
1808
                      end if
1809
                   end do
1810

    
1811
                case (3)
1812
                   !---- Mxyz ----!
1813
                   do j=1,3
1814
                      if (FmAtom%atom(na)%lSkR(j,ik) /=0) then
1815
                         nc=FmAtom%atom(na)%lSkR(j,ik)
1816
                         call Delete_RefCodes(nc,FmAtom)
1817
                      end if
1818
                   end do
1819
                   do j=1,3
1820
                      if (FmAtom%atom(na)%lSkI(j,ik) /=0) then
1821
                         nc=FmAtom%atom(na)%lSkI(j,ik)
1822
                         call Delete_RefCodes(nc,FmAtom)
1823
                      end if
1824
                   end do
1825

    
1826
                case (4:)
1827
                   err_refcodes=.true.
1828
                   ERR_RefCodes_Mess="Incompatible information for this type of variable "
1829
                   return
1830
             end select
1831

    
1832
          !---- VARY Directive ----!
1833
          case ("var")
1834

    
1835
              select case (key)
1836
                case (0)
1837

    
1838
                   !---- nb must be different zero ----!
1839
                   select case (nb)
1840
                      case (0)
1841
                         err_refcodes=.true.
1842
                         ERR_RefCodes_Mess="Option not defined"
1843
                         return
1844

    
1845
                      case ( 1:3)
1846
                         !---- Rx_, Ry_, Rz_ ----!
1847
                         if (FmAtom%atom(na)%lSkR(nb,ik) ==0) then
1848

    
1849
                            np_refi=np_refi+1
1850
                            FmAtom%atom(na)%lSkR(nb,ik)=np_refi
1851
                            FmAtom%atom(na)%mSkR(nb,ik)=1.0
1852

    
1853
                            if (FmAtom%atom(na)%lSkR(nb,ik) == np_refi) then
1854
                               V_Vec(np_refi)=FmAtom%atom(na)%SkR(nb,ik)
1855
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1856
                               V_Bounds(1,np_refi)=xl
1857
                               V_Bounds(2,np_refi)=xu
1858
                               V_Bounds(3,np_refi)=xs
1859
                               V_BCon(np_refi)=ic
1860
                               V_list(np_refi)=na
1861
                            else
1862
                               np_refi=np_refi-1
1863
                            end if
1864
                         end if
1865

    
1866
                      case ( 4:6)
1867
                         !---- Ix_, Iy_, Iz_ ----!
1868
                         if (FmAtom%atom(na)%lSkI(nb-3,ik) ==0) then
1869

    
1870
                            np_refi=np_refi+1
1871
                            FmAtom%atom(na)%lSkI(nb-3,ik)=np_refi
1872
                            FmAtom%atom(na)%mSkI(nb-3,ik)=1.0
1873

    
1874
                            if (FmAtom%atom(na)%lSkI(nb-3,ik) == np_refi) then
1875
                               V_Vec(np_refi)=FmAtom%atom(na)%SkI(nb-3,ik)
1876
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1877
                               V_Bounds(1,np_refi)=xl
1878
                               V_Bounds(2,np_refi)=xu
1879
                               V_Bounds(3,np_refi)=xs
1880
                               V_BCon(np_refi)=ic
1881
                               V_list(np_refi)=na
1882
                             else
1883
                               np_refi=np_refi-1
1884
                            end if
1885
                         end if
1886

    
1887
                      case ( 7:9)
1888
                         !---- Rm_, Rphi_, Rth_ ----!
1889
                         if (FmAtom%atom(na)%lSkR(nb-6,ik) ==0) then
1890

    
1891
                            np_refi=np_refi+1
1892
                            FmAtom%atom(na)%lSkR(nb-6,ik)=np_refi
1893
                            FmAtom%atom(na)%mSkR(nb-6,ik)=1.0
1894

    
1895
                            if (FmAtom%atom(na)%lSkR(nb-6,ik) == np_refi) then
1896
                               V_Vec(np_refi)=FmAtom%atom(na)%Spher_SkR(nb-6,ik)
1897
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1898
                               V_Bounds(1,np_refi)=xl
1899
                               V_Bounds(2,np_refi)=xu
1900
                               V_Bounds(3,np_refi)=xs
1901
                               V_BCon(np_refi)=ic
1902
                               V_list(np_refi)=na
1903
                            else
1904
                               np_refi=np_refi-1
1905
                            end if
1906
                         end if
1907

    
1908
                      case (10:12)
1909
                         !---- Im_, Iphi_, Ith_ ----!
1910
                         if (FmAtom%atom(na)%lSkI(nb-9,ik) ==0) then
1911

    
1912
                            np_refi=np_refi+1
1913
                            FmAtom%atom(na)%lSkI(nb-9,ik)=np_refi
1914
                            FmAtom%atom(na)%mSkI(nb-9,ik)=1.0
1915

    
1916
                            if (FmAtom%atom(na)%lSkI(nb-9,ik) == np_refi) then
1917
                               V_Vec(np_refi)=FmAtom%atom(na)%Spher_SkI(nb-9,ik)
1918
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1919
                               V_Bounds(1,np_refi)=xl
1920
                               V_Bounds(2,np_refi)=xu
1921
                               V_Bounds(3,np_refi)=xs
1922
                               V_BCon(np_refi)=ic
1923
                               V_list(np_refi)=na
1924
                             else
1925
                               np_refi=np_refi-1
1926
                            end if
1927
                         end if
1928

    
1929
                      case (13)
1930
                         !---- MagPh_ ----!
1931
                         if (FmAtom%atom(na)%lmphas(ik) ==0) then
1932

    
1933
                            np_refi=np_refi+1
1934
                            FmAtom%atom(na)%mmphas(ik)=1.0
1935
                            FmAtom%atom(na)%lmphas(ik)=np_refi
1936

    
1937
                            V_Vec(np_refi)=FmAtom%atom(na)%mphas(ik)
1938
                            V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1939
                            V_Bounds(1,np_refi)=xl
1940
                            V_Bounds(2,np_refi)=xu
1941
                            V_Bounds(3,np_refi)=xs
1942
                            V_BCon(np_refi)=ic
1943
                            V_list(np_refi)=na
1944
                         end if
1945

    
1946
                      case (14:22)
1947
                         !---- C1_,..C12_ ----!
1948
                         if (FmAtom%atom(na)%lbas(nb-13,ik) ==0) then
1949

    
1950
                            np_refi=np_refi+1
1951
                            FmAtom%atom(na)%lbas(nb-13,ik)=np_refi
1952
                            FmAtom%atom(na)%mbas(nb-13,ik)=1.0
1953

    
1954
                            if (FmAtom%atom(na)%lbas(nb-13,ik) == np_refi) then
1955
                               V_Vec(np_refi)=FmAtom%atom(na)%cbas(nb-13,ik)
1956
                               V_Name(np_refi)=trim(mcode_nam(nb))//trim(FmAtom%atom(na)%lab)
1957
                               V_Bounds(1,np_refi)=xl
1958
                               V_Bounds(2,np_refi)=xu
1959
                               V_Bounds(3,np_refi)=xs
1960
                               V_BCon(np_refi)=ic
1961
                               V_list(np_refi)=na
1962
                             else
1963
                               np_refi=np_refi-1
1964
                            end if
1965
                         end if
1966

    
1967
                      case (23:)
1968
                         err_refcodes=.true.
1969
                         ERR_RefCodes_Mess="Option Not defined by this type of variables"
1970
                         return
1971

    
1972
                   end select ! nb
1973

    
1974
                case (1)
1975
                   !---- Rxyz ----!
1976
                   do j=1,3
1977
                      if (FmAtom%atom(na)%lSkR(j,ik) ==0) then
1978

    
1979
                            np_refi=np_refi+1
1980
                            FmAtom%atom(na)%lSkR(j,ik)=np_refi
1981
                            FmAtom%atom(na)%mSkR(j,ik)=1.0
1982

    
1983
                         if (FmAtom%atom(na)%lSkR(j,ik) == np_refi) then
1984
                            V_Vec(np_refi)=FmAtom%atom(na)%SkR(j,ik)
1985
                            V_Name(np_refi)=trim(mcode_nam(j))//trim(FmAtom%atom(na)%lab)
1986
                            V_Bounds(1,np_refi)=xl
1987
                            V_Bounds(2,np_refi)=xu
1988
                            V_Bounds(3,np_refi)=xs
1989
                            V_BCon(np_refi)=ic
1990
                            V_list(np_refi)=na
1991
                         else
1992
                            np_refi=np_refi-1
1993
                         end if
1994
                      end if
1995
                   end do
1996

    
1997
                case (2)
1998
                   !---- Ixyz ----!
1999
                   do j=1,3
2000
                      if (FmAtom%atom(na)%lSkI(j,ik) ==0) then
2001

    
2002
                            np_refi=np_refi+1
2003
                            FmAtom%atom(na)%lSkI(j,ik)=np_refi
2004
                            FmAtom%atom(na)%mSkI(j,ik)=1.0
2005

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

    
2020
                case (3)
2021
                    !---- Mxyz ----!
2022
                   do j=1,3
2023
                      if (FmAtom%atom(na)%lSkR(j,ik) ==0) then
2024

    
2025
                            np_refi=np_refi+1
2026
                            FmAtom%atom(na)%lSkR(j,ik)=np_refi
2027
                            FmAtom%atom(na)%mSkR(j,ik)=1.0
2028

    
2029
                         if (FmAtom%atom(na)%lSkR(j,ik) == np_refi) then
2030
                            V_Vec(np_refi)=FmAtom%atom(na)%SkR(j,ik)
2031
                            V_Name(np_refi)=trim(mcode_nam(j))//trim(FmAtom%atom(na)%lab)
2032
                            V_Bounds(1,np_refi)=xl
2033
                            V_Bounds(2,np_refi)=xu
2034
                            V_Bounds(3,np_refi)=xs
2035
                            V_BCon(np_refi)=ic
2036
                            V_list(np_refi)=na
2037
                         else
2038
                            np_refi=np_refi-1
2039
                         end if
2040
                      end if
2041
                   end do
2042
                   
2043
                   do j=1,3
2044
                      if (FmAtom%atom(na)%lSkI(j,ik) ==0) then
2045

    
2046
                            np_refi=np_refi+1
2047
                            FmAtom%atom(na)%lSkI(j,ik)=np_refi
2048
                            FmAtom%atom(na)%mSkI(j,ik)=1.0
2049

    
2050
                         if (FmAtom%atom(na)%lSkI(j,ik) == np_refi) then
2051
                            V_Vec(np_refi)=FmAtom%atom(na)%SkI(j,ik)
2052
                            V_Name(np_refi)=trim(mcode_nam(j+3))//trim(FmAtom%atom(na)%lab)
2053
                            V_Bounds(1,np_refi)=xl
2054
                            V_Bounds(2,np_refi)=xu
2055
                            V_Bounds(3,np_refi)=xs
2056
                            V_BCon(np_refi)=ic
2057
                            V_list(np_refi)=na
2058
                         else
2059
                            np_refi=np_refi-1
2060
                         end if
2061
                      end if
2062
                   end do
2063
 
2064
                case (4:)
2065
                   err_refcodes=.true.
2066
                   ERR_RefCodes_Mess="Option Not defined by this type of variables"
2067
                   return
2068
             end select
2069
       end select
2070

    
2071
       return
2072
    End Subroutine Fill_RefCodes_FmAtom
2073

    
2074
    !!--++
2075
    !!--++ Subroutine Fill_RefCodes_Molcrys(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molcrys,NMol)
2076
    !!--++    integer,                      intent(in)     :: Key
2077
    !!--++    character(len=*),             intent(in)     :: Dire
2078
    !!--++    integer,                      intent(in)     :: Na
2079
    !!--++    integer,                      intent(in)     :: Nb
2080
    !!--++    real(kind=cp),                intent(in)     :: Xl
2081
    !!--++    real(kind=cp),                intent(in)     :: Xu
2082
    !!--++    real(kind=cp),                intent(in)     :: Xs
2083
    !!--++    integer,                      intent(in)     :: Ic
2084
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
2085
    !!--++    integer,                      intent(in)     :: NMol
2086
    !!--++
2087
    !!--++ Overloaded
2088
    !!--++ Write on Vectors the Information for Free Atoms
2089
    !!--++
2090
    !!--++ Update: March - 2005
2091
    !!
2092
    Subroutine Fill_RefCodes_Molcrys(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molcrys,Nmol)
2093
       !---- Arguments ----!
2094
       integer,                      intent(in)     :: Key
2095
       character(len=*),             intent(in)     :: Dire
2096
       integer,                      intent(in)     :: Na
2097
       integer,                      intent(in)     :: Nb
2098
       real(kind=cp),                intent(in)     :: Xl
2099
       real(kind=cp),                intent(in)     :: Xu
2100
       real(kind=cp),                intent(in)     :: Xs
2101
       integer,                      intent(in)     :: Ic
2102
       type(molecular_Crystal_type), intent(in out) :: MolCrys
2103
       integer,                      intent(in)     :: NMol
2104

    
2105
       !---- Local variables ----!
2106
       character(len=2) :: car
2107
       integer          :: i, j, nc, naa
2108

    
2109
       call init_err_refcodes()
2110
       if (Na <= 0) then
2111
          err_refcodes=.true.
2112
          ERR_RefCodes_Mess="Number of atom no defined"
2113
          return
2114
       end if
2115

    
2116
       select case (dire)
2117
          !---- FIX Directive ----!
2118
          case ("fix")
2119

    
2120
             select case (key)
2121
                case (0)
2122
                   !---- Nb must be different zero ----!
2123
                   select case (nb)
2124
                      case (0)
2125
                         err_refcodes=.true.
2126
                         ERR_RefCodes_Mess="Option Not defined"
2127
                         return
2128

    
2129
                      case ( 1:3)
2130
                         !---- X_, Y_, Z_ ----!
2131
                         if (na <= molcrys%n_free) then
2132
                            if (molcrys%atm(na)%lx(nb) /=0) then
2133
                               nc=molcrys%atm(na)%lx(nb)
2134
                               call Delete_RefCodes(nc,MolCrys)
2135
                            end if
2136
                         else
2137
                            naa=na-molcrys%n_free
2138
                            do i=1,molcrys%n_mol
2139
                               if (naa > molcrys%mol(i)%natoms) then
2140
                                  naa=naa-molcrys%mol(i)%natoms
2141
                                  cycle
2142
                               end if
2143
                               if (molcrys%mol(i)%lI_coor(nb,naa) /=0) then
2144
                                  nc=molcrys%mol(i)%lI_coor(nb,naa)
2145
                                  call Delete_RefCodes(nc,MolCrys)
2146
                               end if
2147
                            end do
2148
                         end if
2149

    
2150
                      case ( 4)
2151
                         !---- Biso_ ----!
2152
                         if (na <= molcrys%n_free) then
2153
                            if (molcrys%atm(na)%lbiso /=0) then
2154
                               nc=molcrys%atm(na)%lbiso
2155
                               call Delete_RefCodes(nc,MolCrys)
2156
                            end if
2157
                         else
2158
                            naa=na-molcrys%n_free
2159
                            do i=1,molcrys%n_mol
2160
                               if (naa > molcrys%mol(i)%natoms) then
2161
                                  naa=naa-molcrys%mol(i)%natoms
2162
                                  cycle
2163
                               end if
2164
                               if (molcrys%mol(i)%lbiso(naa) /=0) then
2165
                                  nc=molcrys%mol(i)%lbiso(naa)
2166
                                  call Delete_RefCodes(nc,MolCrys)
2167
                               end if
2168
                            end do
2169
                         end if
2170

    
2171
                      case ( 5)
2172
                         !---- Occ_ ----!
2173
                         if (na <= molcrys%n_free) then
2174
                            if (molcrys%atm(na)%locc /=0) then
2175
                               nc=molcrys%atm(na)%locc
2176
                               call Delete_RefCodes(nc,MolCrys)
2177
                            end if
2178
                         else
2179
                            naa=na-molcrys%n_free
2180
                            do i=1,molcrys%n_mol
2181
                               if (naa > molcrys%mol(i)%natoms) then
2182
                                  naa=naa-molcrys%mol(i)%natoms
2183
                                  cycle
2184
                               end if
2185
                               if (molcrys%mol(i)%locc(naa) /=0) then
2186
                                  nc=molcrys%mol(i)%locc(naa)
2187
                                  call Delete_RefCodes(nc,MolCrys)
2188
                               end if
2189
                            end do
2190
                         end if
2191

    
2192
                      case ( 6:11)
2193
                         !---- B11_, ..., B23_ ----!
2194
                         if (na <= molcrys%n_free) then
2195
                            if (molcrys%atm(na)%lu(nb-5) /=0) then
2196
                               nc=molcrys%atm(na)%lu(nb-5)
2197
                               call Delete_RefCodes(nc,MolCrys)
2198
                            end if
2199
                         else
2200
                            err_refcodes=.true.
2201
                            ERR_RefCodes_Mess="Option Not defined"
2202
                            return
2203
                         end if
2204

    
2205
                      case (12)
2206
                         !---- Banis_ ----!
2207
                         if (na <= molcrys%n_free) then
2208
                            do j=1,6
2209
                               if (molcrys%atm(na)%lu(j) /=0) then
2210
                                  nc=molcrys%atm(na)%lu(j)
2211
                                  call Delete_RefCodes(nc,MolCrys)
2212
                               end if
2213
                            end do
2214
                         else
2215
                            err_refcodes=.true.
2216
                            ERR_RefCodes_Mess="Option Not defined"
2217
                            return
2218
                         end if
2219

    
2220
                      case (13:15)
2221
                         !---- Xc_, Yc_, Zc_ ----!
2222
                         select case (nmol)
2223
                            case (-1)
2224
                               err_refcodes=.true.
2225
                               ERR_RefCodes_Mess="Option Not defined"
2226
                               return
2227

    
2228
                            case (0)
2229
                               do i=1,molcrys%n_mol
2230
                                  if (molcrys%mol(i)%lxcentre(nb-12) /=0) then
2231
                                     nc=molcrys%mol(i)%lxcentre(nb-12)
2232
                                     call Delete_RefCodes(nc,MolCrys)
2233
                                  end if
2234
                               end do
2235

    
2236
                            case (1:)
2237
                               if (molcrys%mol(nmol)%lxcentre(nb-12) /=0) then
2238
                                  nc=molcrys%mol(nmol)%lxcentre(nb-12)
2239
                                  call Delete_RefCodes(nc,MolCrys)
2240
                               end if
2241
                         end select
2242

    
2243
                      case (16:18)
2244
                         !---- Theta_, Phi_, Chi_ ----!
2245
                         select case (nmol)
2246
                            case (-1)
2247
                               err_refcodes=.true.
2248
                               ERR_RefCodes_Mess="Option Not defined"
2249
                               return
2250

    
2251
                            case (0)
2252
                               do i=1,molcrys%n_mol
2253
                                  if (molcrys%mol(i)%lOrient(nb-15) /=0) then
2254
                                     nc=molcrys%mol(i)%lOrient(nb-15)
2255
                                     call Delete_RefCodes(nc,MolCrys)
2256
                                  end if
2257
                               end do
2258

    
2259
                            case (1:)
2260
                               if (molcrys%mol(nmol)%lOrient(nb-15) /=0) then
2261
                                  nc=molcrys%mol(nmol)%lOrient(nb-15)
2262
                                  call Delete_RefCodes(nc,MolCrys)
2263
                               end if
2264
                         end select
2265

    
2266
                      case (19:21)
2267
                         !!! Not yet implemented !!!
2268

    
2269
                   end select ! nb
2270

    
2271
                case (1)
2272
                   !---- XYZ ----!
2273
                   if (na <= molcrys%n_free) then
2274
                      do j=1,3
2275
                         if (molcrys%atm(na)%lx(j) /=0) then
2276
                            nc=molcrys%atm(na)%lx(j)
2277
                            call Delete_RefCodes(nc,MolCrys)
2278
                         end if
2279
                      end do
2280
                   else
2281
                      naa=na-molcrys%n_free
2282
                      do i=1,molcrys%n_mol
2283
                         if (naa > molcrys%mol(i)%natoms) then
2284
                            naa=naa-molcrys%mol(i)%natoms
2285
                            cycle
2286
                         end if
2287
                         do j=1,3
2288
                            if (molcrys%mol(i)%lI_coor(j,naa) /=0) then
2289
                               nc=molcrys%mol(i)%lI_coor(j,naa)
2290
                               call Delete_RefCodes(nc,MolCrys)
2291
                            end if
2292
                         end do
2293
                      end do
2294
                   end if
2295

    
2296
                case (2)
2297
                   !---- OCC ----!
2298
                   if (na <= molcrys%n_free) then
2299
                      if (molcrys%atm(na)%locc /=0) then
2300
                         nc=molcrys%atm(na)%locc
2301
                         call Delete_RefCodes(nc,MolCrys)
2302
                      end if
2303
                   else
2304
                      naa=na-molcrys%n_free
2305
                      do i=1,molcrys%n_mol
2306
                         if (naa > molcrys%mol(i)%natoms) then
2307
                            naa=naa-molcrys%mol(i)%natoms
2308
                            cycle
2309
                         end if
2310
                         if (molcrys%mol(i)%locc(naa) /=0) then
2311
                            nc=molcrys%mol(i)%locc(naa)
2312
                            call Delete_RefCodes(nc,MolCrys)
2313
                         end if
2314
                      end do
2315
                   end if
2316

    
2317
                case (3)
2318
                   !---- BIS ----!
2319
                   if (na <= molcrys%n_free) then
2320
                      if (molcrys%atm(na)%lbiso /=0) then
2321
                         nc=molcrys%atm(na)%lbiso
2322
                         call Delete_RefCodes(nc,MolCrys)
2323
                      end if
2324
                   else
2325
                      naa=na-molcrys%n_free
2326
                      do i=1,molcrys%n_mol
2327
                         if (naa > molcrys%mol(i)%natoms) then
2328
                            naa=naa-molcrys%mol(i)%natoms
2329
                            cycle
2330
                         end if
2331
                         if (molcrys%mol(i)%lbiso(naa) /=0) then
2332
                            nc=molcrys%mol(i)%lbiso(naa)
2333
                            call Delete_RefCodes(nc,MolCrys)
2334
                         end if
2335
                      end do
2336
                   end if
2337

    
2338
                case (4)
2339
                   !---- BAN ----!
2340
                   if (na <= molcrys%n_free) then
2341
                      do j=1,6
2342
                         if (molcrys%atm(na)%lu(j) /=0) then
2343
                            nc=molcrys%atm(na)%lu(j)
2344
                            call Delete_RefCodes(nc,MolCrys)
2345
                         end if
2346
                      end do
2347
                   else
2348
                      err_refcodes=.true.
2349
                      ERR_RefCodes_Mess="Option Not defined"
2350
                      return
2351
                   end if
2352

    
2353
                case (5)
2354
                   !---- ALL ----!
2355
                   if (na <= molcrys%n_free) then
2356
                      do j=1,3
2357
                         if (molcrys%atm(na)%lx(j) /=0) then
2358
                            nc=molcrys%atm(na)%lx(j)
2359
                            call Delete_RefCodes(nc,MolCrys)
2360
                         end if
2361
                      end do
2362
                   else
2363
                      naa=na-molcrys%n_free
2364
                      do i=1,molcrys%n_mol
2365
                         if (naa > molcrys%mol(i)%natoms) then
2366
                            naa=naa-molcrys%mol(i)%natoms
2367
                            cycle
2368
                         end if
2369
                         do j=1,3
2370
                            if (molcrys%mol(i)%lI_coor(j,naa) /=0) then
2371
                               nc=molcrys%mol(i)%lI_coor(j,naa)
2372
                               call Delete_RefCodes(nc,MolCrys)
2373
                            end if
2374
                         end do
2375
                      end do
2376
                   end if
2377

    
2378
                   if (na <= molcrys%n_free) then
2379
                      if (molcrys%atm(na)%locc /=0) then
2380
                         nc=molcrys%atm(na)%locc
2381
                         call Delete_RefCodes(nc,MolCrys)
2382
                      end if
2383
                   else
2384
                      naa=na-molcrys%n_free
2385
                      do i=1,molcrys%n_mol
2386
                         if (naa > molcrys%mol(i)%natoms) then
2387
                            naa=naa-molcrys%mol(i)%natoms
2388
                            cycle
2389
                         end if
2390
                         if (molcrys%mol(i)%locc(naa) /=0) then
2391
                            nc=molcrys%mol(i)%locc(naa)
2392
                            call Delete_RefCodes(nc,MolCrys)
2393
                         end if
2394
                      end do
2395
                   end if
2396

    
2397
                   if (na <= molcrys%n_free) then
2398
                      if (molcrys%atm(na)%lbiso /=0) then
2399
                         nc=molcrys%atm(na)%lbiso
2400
                         call Delete_RefCodes(nc,MolCrys)
2401
                      end if
2402
                   else
2403
                      naa=na-molcrys%n_free
2404
                      do i=1,molcrys%n_mol
2405
                         if (naa > molcrys%mol(i)%natoms) then
2406
                            naa=naa-molcrys%mol(i)%natoms
2407
                            cycle
2408
                         end if
2409
                         if (molcrys%mol(i)%lbiso(naa) /=0) then
2410
                            nc=molcrys%mol(i)%lbiso(naa)
2411
                            call Delete_RefCodes(nc,MolCrys)
2412
                         end if
2413
                      end do
2414
                   end if
2415

    
2416
                   if (na <= molcrys%n_free) then
2417
                      do j=1,6
2418
                         if (molcrys%atm(na)%lu(j) /=0) then
2419
                            nc=molcrys%atm(na)%lu(j)
2420
                            call Delete_RefCodes(nc,MolCrys)
2421
                         end if
2422
                      end do
2423
                   end if
2424

    
2425
                   select case (nmol)
2426
                      case (-1)
2427
                         err_refcodes=.true.
2428
                         ERR_RefCodes_Mess="Option Not defined"
2429
                         return
2430

    
2431
                      case (0)
2432
                         do i=1,molcrys%n_mol
2433
                            do j=1,3
2434
                               if (molcrys%mol(i)%lxcentre(j) /=0) then
2435
                                  nc=molcrys%mol(i)%lxcentre(j)
2436
                                  call Delete_RefCodes(nc,molcrys)
2437
                               end if
2438
                            end do
2439

    
2440
                            do j=1,3
2441
                               if (molcrys%mol(i)%lOrient(j) /=0) then
2442
                                  nc=molcrys%mol(i)%lOrient(j)
2443
                                  call Delete_RefCodes(nc,molcrys)
2444
                               end if
2445
                            end do
2446
                         end do
2447

    
2448
                      case (1:)
2449
                         do j=1,3
2450
                            if (molcrys%mol(nmol)%lxcentre(j) /=0) then
2451
                               nc=molcrys%mol(nmol)%lxcentre(j)
2452
                               call Delete_RefCodes(nc,molcrys)
2453
                            end if
2454
                         end do
2455

    
2456
                         do j=1,3
2457
                            if (molcrys%mol(nmol)%lOrient(j) /=0) then
2458
                               nc=molcrys%mol(nmol)%lOrient(j)
2459
                               call Delete_RefCodes(nc,molcrys)
2460
                            end if
2461
                         end do
2462
                   end select
2463

    
2464
                   !!! Not yet Implemented !!!
2465

    
2466
                case (6)
2467
                   !---- CEN ----!
2468
                   select case (nmol)
2469
                      case (-1)
2470
                         err_refcodes=.true.
2471
                         ERR_RefCodes_Mess="Option Not defined"
2472
                         return
2473

    
2474
                      case (0)
2475
                         do i=1,molcrys%n_mol
2476
                            do j=1,3
2477
                               if (molcrys%mol(i)%lxcentre(j) /=0) then
2478
                                  nc=molcrys%mol(i)%lxcentre(j)
2479
                                  call Delete_RefCodes(nc,molcrys)
2480
                               end if
2481
                            end do
2482
                         end do
2483

    
2484
                      case (1:)
2485
                         do j=1,3
2486
                            if (molcrys%mol(nmol)%lxcentre(j) /=0) then
2487
                               nc=molcrys%mol(nmol)%lxcentre(j)
2488
                               call Delete_RefCodes(nc,molcrys)
2489
                            end if
2490
                         end do
2491
                   end select
2492

    
2493
                case (7)
2494
                   !---- ORI ----!
2495
                   select case (nmol)
2496
                      case (-1)
2497
                         err_refcodes=.true.
2498
                         ERR_RefCodes_Mess="Option Not defined"
2499
                         return
2500

    
2501
                      case (0)
2502
                         do i=1,molcrys%n_mol
2503
                            do j=1,3
2504
                               if (molcrys%mol(i)%lOrient(j) /=0) then
2505
                                  nc=molcrys%mol(i)%lOrient(j)
2506
                                  call Delete_RefCodes(nc,molcrys)
2507
                               end if
2508
                            end do
2509
                         end do
2510

    
2511
                      case (1:)
2512
                         do j=1,3
2513
                            if (molcrys%mol(nmol)%lOrient(j) /=0) then
2514
                              nc=molcrys%mol(nmol)%lOrient(j)
2515
                               call Delete_RefCodes(nc,molcrys)
2516
                            end if
2517
                         end do
2518
                   end select
2519

    
2520
                case (8)
2521
                   !---- THE ----!
2522
                   !!! Not yet Implemented !!!!
2523
             end select
2524

    
2525
          !---- VARY Directive ----!
2526
          case ("var")
2527

    
2528
             select case (key)
2529
                case (0)
2530
                   !---- Nb must be different zero ----!
2531
                   select case (nb)
2532
                      case (0)
2533
                         err_refcodes=.true.
2534
                         ERR_RefCodes_Mess="Option Not defined"
2535
                         return
2536

    
2537
                      case ( 1:3)
2538
                         !---- X_, Y_, Z_ ----!
2539
                         if (na <= molcrys%n_free) then
2540
                            if (molcrys%atm(na)%lx(nb) ==0) then
2541
                               molcrys%atm(na)%mx(nb)=1.0
2542
                               call get_atompos_ctr(molcrys%atm(na)%x,   &
2543
                                                    molcrys%Spg,np_refi, &
2544
                                                    molcrys%atm(na)%lx,  &
2545
                                                    molcrys%atm(na)%mx)
2546
                               if (molcrys%atm(na)%lx(nb) == np_refi) then
2547
                                  V_Vec(np_refi)=molcrys%atm(na)%x(nb)
2548
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2549
                                  V_Bounds(1,np_refi)=xl
2550
                                  V_Bounds(2,np_refi)=xu
2551
                                  V_Bounds(3,np_refi)=xs
2552
                                  V_BCon(np_refi)=ic
2553
                                  V_list(np_refi)=na
2554
                               else
2555
                                  np_refi=np_refi-1
2556
                               end if
2557
                            end if
2558
                         else
2559
                            naa=na-molcrys%n_free
2560
                            do i=1,molcrys%n_mol
2561
                               if (naa > molcrys%mol(i)%natoms) then
2562
                                  naa=naa-molcrys%mol(i)%natoms
2563
                                  cycle
2564
                               end if
2565

    
2566
                               if (molcrys%mol(i)%lI_coor(nb,naa) ==0) then
2567
                                  molcrys%mol(i)%mI_coor(nb,naa)=1.0
2568
                                  call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa),  &
2569
                                                       molcrys%Spg, np_refi,          &
2570
                                                       molcrys%mol(i)%lI_coor(:,naa), &
2571
                                                       molcrys%mol(i)%mI_coor(:,naa))
2572
                                  if (molcrys%mol(i)%lI_coor(nb,naa) == np_refi) then
2573
                                     V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
2574
                                     V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2575
                                     V_Bounds(1,np_refi)=xl
2576
                                     V_Bounds(2,np_refi)=xu
2577
                                     V_Bounds(3,np_refi)=xs
2578
                                     V_BCon(np_refi)=ic
2579
                                     V_list(np_refi)=na
2580
                                  else
2581
                                     np_refi=np_refi-1
2582
                                  end if
2583
                               end if
2584
                            end do
2585
                         end if
2586

    
2587
                      case ( 4)
2588
                         !---- Biso_ ----!
2589
                         if (na <= molcrys%n_free) then
2590
                            if (molcrys%atm(na)%lbiso ==0) then
2591
                               np_refi=np_refi+1
2592
                               V_Vec(np_refi)=molcrys%atm(na)%biso
2593
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2594
                               molcrys%atm(na)%mbiso=1.0
2595
                               molcrys%atm(na)%lbiso=np_refi
2596
                               V_Bounds(1,np_refi)=xl
2597
                               V_Bounds(2,np_refi)=xu
2598
                               V_Bounds(3,np_refi)=xs
2599
                               V_BCon(np_refi)=ic
2600
                               V_list(np_refi)=na
2601
                            end if
2602
                         else
2603
                            naa=na-molcrys%n_free
2604
                            do i=1,molcrys%n_mol
2605
                               if (naa > molcrys%mol(i)%natoms) then
2606
                                  naa=naa-molcrys%mol(i)%natoms
2607
                                  cycle
2608
                               end if
2609
                               if (molcrys%mol(i)%lbiso(naa) ==0) then
2610
                                  np_refi=np_refi+1
2611
                                  V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
2612
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2613
                                  molcrys%mol(i)%mbiso(naa)=1.0
2614
                                  molcrys%mol(i)%lbiso(naa)=np_refi
2615
                                  V_Bounds(1,np_refi)=xl
2616
                                  V_Bounds(2,np_refi)=xu
2617
                                  V_Bounds(3,np_refi)=xs
2618
                                  V_BCon(np_refi)=ic
2619
                                  V_list(np_refi)=na
2620
                               end if
2621
                            end do
2622
                         end if
2623

    
2624
                      case ( 5)
2625
                         !---- Occ_ ----!
2626
                         if (na <= molcrys%n_free) then
2627
                            if (molcrys%atm(na)%locc ==0) then
2628
                               np_refi=np_refi+1
2629
                               V_Vec(np_refi)=molcrys%atm(na)%occ
2630
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2631
                               molcrys%atm(na)%mocc=1.0
2632
                               molcrys%atm(na)%locc=np_refi
2633
                               V_Bounds(1,np_refi)=xl
2634
                               V_Bounds(2,np_refi)=xu
2635
                               V_Bounds(3,np_refi)=xs
2636
                               V_BCon(np_refi)=ic
2637
                               V_list(np_refi)=na
2638
                            end if
2639
                         else
2640
                            naa=na-molcrys%n_free
2641
                            do i=1,molcrys%n_mol
2642
                               if (naa > molcrys%mol(i)%natoms) then
2643
                                  naa=naa-molcrys%mol(i)%natoms
2644
                                  cycle
2645
                               end if
2646
                               if (molcrys%mol(i)%locc(naa) ==0) then
2647
                                  np_refi=np_refi+1
2648
                                  V_Vec(np_refi)=molcrys%mol(i)%occ(naa)
2649
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%mol(i)%AtName(naa))
2650
                                  molcrys%mol(i)%mocc(naa)=1.0
2651
                                  molcrys%mol(i)%locc(naa)=np_refi
2652
                                  V_Bounds(1,np_refi)=xl
2653
                                  V_Bounds(2,np_refi)=xu
2654
                                  V_Bounds(3,np_refi)=xs
2655
                                  V_BCon(np_refi)=ic
2656
                                  V_list(np_refi)=na
2657
                               end if
2658
                            end do
2659
                         end if
2660

    
2661
                      case ( 6:11)
2662
                         !---- B11_, ..., B23_ ----!
2663
                         if (na <= molcrys%n_free) then
2664
                            if (molcrys%atm(na)%lu(nb-5) ==0) then
2665
                               molcrys%atm(na)%mu(nb-5)=1.0
2666
                               call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
2667
                                                   np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
2668
                               if (molcrys%atm(na)%lu(nb-5) == np_refi) then
2669
                                  V_Vec(np_refi)=molcrys%atm(na)%u(nb-5)
2670
                                  V_Name(np_refi)=trim(code_nam(nb))//trim(molcrys%atm(na)%lab)
2671
                                  V_Bounds(1,np_refi)=xl
2672
                                  V_Bounds(2,np_refi)=xu
2673
                                  V_Bounds(3,np_refi)=xs
2674
                                  V_BCon(np_refi)=ic
2675
                                  V_list(np_refi)=na
2676
                               else
2677
                                  np_refi=np_refi-1
2678
                               end if
2679
                            end if
2680
                         else
2681
                            err_refcodes=.true.
2682
                            ERR_RefCodes_Mess="Option Not defined"
2683
                            return
2684
                         end if
2685

    
2686
                      case (12)
2687
                         !---- Banis_ ----!
2688
                         if (na <= molcrys%n_free) then
2689
                            do j=1,6
2690
                               if (molcrys%atm(na)%lu(j) ==0) then
2691
                                  molcrys%atm(na)%mu(j)=1.0
2692
                                  call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
2693
                                                       np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
2694
                                  if (molcrys%atm(na)%lu(j) == np_refi) then
2695
                                     V_Vec(np_refi)=molcrys%atm(na)%u(j)
2696
                                     V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
2697
                                     V_Bounds(1,np_refi)=xl
2698
                                     V_Bounds(2,np_refi)=xu
2699
                                     V_Bounds(3,np_refi)=xs
2700
                                     V_BCon(np_refi)=ic
2701
                                     V_list(np_refi)=na
2702
                                  else
2703
                                     np_refi=np_refi-1
2704
                                  end if
2705
                               end if
2706
                            end do
2707
                         else
2708
                            err_refcodes=.true.
2709
                            ERR_RefCodes_Mess="Option Not defined"
2710
                            return
2711
                         end if
2712

    
2713
                      case (13:15)
2714
                         !---- Xc_, Yc_, Zc_ ----!
2715
                         select case (nmol)
2716
                            case (-1)
2717
                               err_refcodes=.true.
2718
                               ERR_RefCodes_Mess="Option Not defined"
2719
                               return
2720

    
2721
                            case (0)
2722
                               do i=1,molcrys%n_mol
2723
                                  write(unit=car,fmt="(i2)") i
2724
                                  car=adjustl(car)
2725
                                  if (molcrys%mol(i)%lxcentre(nb-12) ==0) then
2726
                                     np_refi=np_refi+1
2727
                                     V_Vec(np_refi)=molcrys%mol(i)%xcentre(nb-12)
2728
                                     V_Name(np_refi)=trim(code_nam(nb))//"Mol"//trim(car)
2729
                                     molcrys%mol(i)%mxcentre(nb-12)=1.0
2730
                                     molcrys%mol(i)%lxcentre(nb-12)=np_refi
2731
                                     V_Bounds(1,np_refi)=xl
2732
                                     V_Bounds(2,np_refi)=xu
2733
                                     V_Bounds(3,np_refi)=xs
2734
                                     V_BCon(np_refi)=ic
2735
                                     V_list(np_refi)=-i
2736
                                  end if
2737
                               end do
2738

    
2739
                            case (1:)
2740
                               write(unit=car,fmt="(i2)") nmol
2741
                               car=adjustl(car)
2742
                               if (molcrys%mol(nmol)%lxcentre(nb-12) ==0) then
2743
                                  np_refi=np_refi+1
2744
                                  V_Vec(np_refi)=molcrys%mol(nmol)%xcentre(nb-12)
2745
                                  V_Name(np_refi)=trim(code_nam(nb))//"entre_Mol"//trim(car)
2746
                                  molcrys%mol(nmol)%mxcentre(nb-12)=1.0
2747
                                  molcrys%mol(nmol)%lxcentre(nb-12)=np_refi
2748
                                  V_Bounds(1,np_refi)=xl
2749
                                  V_Bounds(2,np_refi)=xu
2750
                                  V_Bounds(3,np_refi)=xs
2751
                                  V_BCon(np_refi)=ic
2752
                                  V_list(np_refi)=-nmol
2753
                               end if
2754
                         end select
2755

    
2756
                      case (16:18)
2757
                         !---- Theta_, Phi_, Chi_ ----!
2758
                         select case (nmol)
2759
                            case (-1)
2760
                               err_refcodes=.true.
2761
                               ERR_RefCodes_Mess="Option Not defined"
2762
                               return
2763

    
2764
                            case (0)
2765
                               do i=1,molcrys%n_mol
2766
                                  write(unit=car,fmt="(i2)") i
2767
                                  car=adjustl(car)
2768
                                  if (molcrys%mol(i)%lOrient(nb-15) ==0) then
2769
                                     np_refi=np_refi+1
2770
                                     V_Vec(np_refi)=molcrys%mol(i)%Orient(nb-15)
2771
                                     V_Name(np_refi)=trim(code_nam(nb))//"Orient_Mol"//trim(car)
2772
                                     molcrys%mol(i)%mOrient(nb-15)=1.0
2773
                                     molcrys%mol(i)%lOrient(nb-15)=np_refi
2774
                                     V_Bounds(1,np_refi)=xl
2775
                                     V_Bounds(2,np_refi)=xu
2776
                                     V_Bounds(3,np_refi)=xs
2777
                                     V_BCon(np_refi)=ic
2778
                                     V_list(np_refi)=-i
2779
                                  end if
2780
                               end do
2781

    
2782
                            case (1:)
2783
                               write(unit=car,fmt="(i2)") nmol
2784
                               car=adjustl(car)
2785
                               if (molcrys%mol(nmol)%lOrient(nb-15) ==0) then
2786
                                  np_refi=np_refi+1
2787
                                  V_Vec(np_refi)=molcrys%mol(nmol)%Orient(nb-15)
2788
                                  V_Name(np_refi)=trim(code_nam(nb))//"Orient_Mol"//trim(car)
2789
                                  molcrys%mol(nmol)%mOrient(nb-15)=1.0
2790
                                  molcrys%mol(nmol)%lOrient(nb-15)=np_refi
2791
                                  V_Bounds(1,np_refi)=xl
2792
                                  V_Bounds(2,np_refi)=xu
2793
                                  V_Bounds(3,np_refi)=xs
2794
                                  V_BCon(np_refi)=ic
2795
                                  V_list(np_refi)=-nmol
2796
                               end if
2797
                         end select
2798

    
2799
                      case (19:21)
2800
                         !!! Not Yet Implemented !!!
2801

    
2802
                   end select ! nb
2803

    
2804
                case (1)
2805
                   !---- XYZ ----!
2806
                   if (na <= molcrys%n_free) then
2807
                      do j=1,3
2808
                         if (molcrys%atm(na)%lx(j) ==0) then
2809
                            molcrys%atm(na)%mx(j)=1.0
2810
                            call get_atompos_ctr(molcrys%atm(na)%x,   &
2811
                                                 molcrys%Spg,np_refi, &
2812
                                                 molcrys%atm(na)%lx,  &
2813
                                                 molcrys%atm(na)%mx)
2814
                            if (molcrys%atm(na)%lx(j) == np_refi) then
2815
                               V_Vec(np_refi)=molcrys%atm(na)%x(j)
2816
                               V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%atm(na)%lab)
2817
                               V_Bounds(1,np_refi)=xl
2818
                               V_Bounds(2,np_refi)=xu
2819
                               V_Bounds(3,np_refi)=xs
2820
                               V_BCon(np_refi)=ic
2821
                               V_list(np_refi)=na
2822
                            else
2823
                               np_refi=np_refi-1
2824
                            end if
2825
                         end if
2826
                      end do
2827
                   else
2828
                      naa=na-molcrys%n_free
2829
                      do i=1,molcrys%n_mol
2830
                         if (naa > molcrys%mol(i)%natoms) then
2831
                             naa=naa-molcrys%mol(i)%natoms
2832
                             cycle
2833
                         end if
2834
                         do j=1,3
2835
                            if (molcrys%mol(i)%lI_coor(j,naa) ==0) then
2836
                               molcrys%mol(i)%mI_coor(nb,naa)=1.0
2837
                               call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa), &
2838
                                                    molcrys%Spg,np_refi,   &
2839
                                                    molcrys%mol(i)%lI_coor(:,naa), &
2840
                                                    molcrys%mol(i)%mI_coor(:,naa))
2841
                               if (molcrys%mol(i)%lI_coor(j,naa) == np_refi) then
2842
                                  V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
2843
                                  V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%mol(i)%AtName(naa))
2844
                                  V_Bounds(1,np_refi)=xl
2845
                                  V_Bounds(2,np_refi)=xu
2846
                                  V_Bounds(3,np_refi)=xs
2847
                                  V_BCon(np_refi)=ic
2848
                                  V_list(np_refi)=na
2849
                               else
2850
                                  np_refi=np_refi-1
2851
                               end if
2852
                            end if
2853
                         end do
2854
                      end do
2855
                   end if
2856

    
2857
                case (2)
2858
                   !---- OCC ----!
2859
                   if (na <= molcrys%n_free) then
2860
                      if (molcrys%atm(na)%locc ==0) then
2861
                         np_refi=np_refi+1
2862
                         V_Vec(np_refi)=molcrys%atm(na)%occ
2863
                         V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%atm(na)%lab)
2864
                         molcrys%atm(na)%mocc=1.0
2865
                         molcrys%atm(na)%locc=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)=na
2871
                      end if
2872
                   else
2873
                      naa=na-molcrys%n_free
2874
                      do i=1,molcrys%n_mol
2875
                         if (naa > molcrys%mol(i)%natoms) then
2876
                            naa=naa-molcrys%mol(i)%natoms
2877
                            cycle
2878
                         end if
2879

    
2880
                         if (molcrys%mol(i)%locc(naa) ==0) then
2881
                            np_refi=np_refi+1
2882
                            V_Vec(np_refi)=molcrys%mol(i)%Occ(naa)
2883
                            V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%mol(i)%AtName(naa))
2884
                            molcrys%mol(i)%mocc(naa)=1.0
2885
                            molcrys%mol(i)%locc(naa)=np_refi
2886
                            V_Bounds(1,np_refi)=xl
2887
                            V_Bounds(2,np_refi)=xu
2888
                            V_Bounds(3,np_refi)=xs
2889
                            V_BCon(np_refi)=ic
2890
                            V_list(np_refi)=na
2891
                         end if
2892
                      end do
2893
                   end if
2894

    
2895
                case (3)
2896
                   !---- BIS ----!
2897
                   if (na <= molcrys%n_free) then
2898
                      if (molcrys%atm(na)%lbiso ==0) then
2899
                         np_refi=np_refi+1
2900
                         V_Vec(np_refi)=molcrys%atm(na)%biso
2901
                         V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%atm(na)%lab)
2902
                         molcrys%atm(na)%mbiso=1.0
2903
                         molcrys%atm(na)%lbiso=np_refi
2904
                         V_Bounds(1,np_refi)=xl
2905
                         V_Bounds(2,np_refi)=xu
2906
                         V_Bounds(3,np_refi)=xs
2907
                         V_BCon(np_refi)=ic
2908
                         V_list(np_refi)=na
2909
                      end if
2910
                   else
2911
                      naa=na-molcrys%n_free
2912
                      do i=1,molcrys%n_mol
2913
                         if (naa > molcrys%mol(i)%natoms) then
2914
                            naa=naa-molcrys%mol(i)%natoms
2915
                            cycle
2916
                         end if
2917
                         if (molcrys%mol(i)%lbiso(naa) ==0) then
2918
                            np_refi=np_refi+1
2919
                            V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
2920
                            V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%mol(i)%AtName(naa))
2921
                            molcrys%mol(i)%mbiso(naa)=1.0
2922
                            molcrys%mol(i)%lbiso(naa)=np_refi
2923
                            V_Bounds(1,np_refi)=xl
2924
                            V_Bounds(2,np_refi)=xu
2925
                            V_Bounds(3,np_refi)=xs
2926
                            V_BCon(np_refi)=ic
2927
                            V_list(np_refi)=na
2928
                         end if
2929
                      end do
2930
                   end if
2931

    
2932
                case (4)
2933
                   !---- BAN ----!
2934
                   if (na <= molcrys%n_free) then
2935
                      do j=1,6
2936
                         if (molcrys%atm(na)%lu(j) ==0) then
2937
                            molcrys%atm(na)%mu(j)=1.0
2938
                            call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
2939
                                                 np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
2940
                            if (molcrys%atm(na)%lu(j) == np_refi) then
2941
                               V_Vec(np_refi)=molcrys%atm(na)%u(j)
2942
                               V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
2943
                               V_Bounds(1,np_refi)=xl
2944
                               V_Bounds(2,np_refi)=xu
2945
                               V_Bounds(3,np_refi)=xs
2946
                               V_BCon(np_refi)=ic
2947
                               V_list(np_refi)=na
2948
                            else
2949
                               np_refi=np_refi-1
2950
                            end if
2951
                         end if
2952
                      end do
2953
                   else
2954
                      err_refcodes=.true.
2955
                      ERR_RefCodes_Mess="Option Not defined"
2956
                      return
2957
                   end if
2958

    
2959
                case (5)
2960
                   !---- ALL ----!
2961
                   if (na <= molcrys%n_free) then
2962
                      do j=1,3
2963
                         if (molcrys%atm(na)%lx(j) ==0) then
2964
                            molcrys%atm(na)%mx(j)=1.0
2965
                            call get_atompos_ctr(molcrys%atm(na)%x,   &
2966
                                                 molcrys%Spg,np_refi, &
2967
                                                 molcrys%atm(na)%lx,  &
2968
                                                 molcrys%atm(na)%mx)
2969
                            if (molcrys%atm(na)%lx(j) == np_refi) then
2970
                               V_Vec(np_refi)=molcrys%atm(na)%x(j)
2971
                               V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%atm(na)%lab)
2972
                               V_Bounds(1,np_refi)=xl
2973
                               V_Bounds(2,np_refi)=xu
2974
                               V_Bounds(3,np_refi)=xs
2975
                               V_BCon(np_refi)=ic
2976
                               V_list(np_refi)=na
2977
                            else
2978
                               np_refi=np_refi-1
2979
                            end if
2980
                         end if
2981
                      end do
2982
                   else
2983
                      naa=na-molcrys%n_free
2984
                      do i=1,molcrys%n_mol
2985
                         if (naa > molcrys%mol(i)%natoms) then
2986
                             naa=naa-molcrys%mol(i)%natoms
2987
                             cycle
2988
                         end if
2989

    
2990
                         do j=1,3
2991
                            if (molcrys%mol(i)%lI_coor(j,naa) ==0) then
2992
                               molcrys%mol(i)%mI_coor(nb,naa)=1.0
2993
                               call get_atompos_ctr(molcrys%mol(i)%I_Coor(:,naa), &
2994
                                                    molcrys%Spg,np_refi,   &
2995
                                                    molcrys%mol(i)%lI_coor(:,naa), &
2996
                                                    molcrys%mol(i)%mI_coor(:,naa))
2997
                               if (molcrys%mol(i)%lI_coor(j,naa) == np_refi) then
2998
                                  V_Vec(np_refi)=molcrys%mol(i)%I_Coor(nb,naa)
2999
                                  V_Name(np_refi)=trim(code_nam(j))//trim(molcrys%mol(i)%AtName(naa))
3000
                                  V_Bounds(1,np_refi)=xl
3001
                                  V_Bounds(2,np_refi)=xu
3002
                                  V_Bounds(3,np_refi)=xs
3003
                                  V_BCon(np_refi)=ic
3004
                                  V_list(np_refi)=na
3005
                               else
3006
                                  np_refi=np_refi-1
3007
                               end if
3008
                            end if
3009
                         end do
3010
                      end do
3011
                   end if
3012

    
3013
                   if (na <= molcrys%n_free) then
3014
                      if (molcrys%atm(na)%locc ==0) then
3015
                         np_refi=np_refi+1
3016
                         V_Vec(np_refi)=molcrys%atm(na)%occ
3017
                         V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%atm(na)%lab)
3018
                         molcrys%atm(na)%mocc=1.0
3019
                         molcrys%atm(na)%locc=np_refi
3020
                         V_Bounds(1,np_refi)=xl
3021
                         V_Bounds(2,np_refi)=xu
3022
                         V_Bounds(3,np_refi)=xs
3023
                         V_BCon(np_refi)=ic
3024
                         V_list(np_refi)=na
3025
                      end if
3026
                   else
3027
                      naa=na-molcrys%n_free
3028
                      do i=1,molcrys%n_mol
3029
                         if (naa > molcrys%mol(i)%natoms) then
3030
                            naa=naa-molcrys%mol(i)%natoms
3031
                            cycle
3032
                         end if
3033
                         if (molcrys%mol(i)%locc(naa) ==0) then
3034
                            np_refi=np_refi+1
3035
                            V_Vec(np_refi)=molcrys%mol(i)%Occ(naa)
3036
                            V_Name(np_refi)=trim(code_nam(5))//trim(molcrys%mol(i)%AtName(naa))
3037
                            molcrys%mol(i)%mocc(naa)=1.0
3038
                            molcrys%mol(i)%locc(naa)=np_refi
3039
                            V_Bounds(1,np_refi)=xl
3040
                            V_Bounds(2,np_refi)=xu
3041
                            V_Bounds(3,np_refi)=xs
3042
                            V_BCon(np_refi)=ic
3043
                            V_list(np_refi)=na
3044
                         end if
3045
                      end do
3046
                   end if
3047

    
3048
                   if (na <= molcrys%n_free) then
3049
                      if (molcrys%atm(na)%lbiso ==0) then
3050
                         np_refi=np_refi+1
3051
                         V_Vec(np_refi)=molcrys%atm(na)%biso
3052
                         V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%atm(na)%lab)
3053
                         molcrys%atm(na)%mbiso=1.0
3054
                         molcrys%atm(na)%lbiso=np_refi
3055
                         V_Bounds(1,np_refi)=xl
3056
                         V_Bounds(2,np_refi)=xu
3057
                         V_Bounds(3,np_refi)=xs
3058
                         V_BCon(np_refi)=ic
3059
                         V_list(np_refi)=na
3060
                      end if
3061
                   else
3062
                      naa=na-molcrys%n_free
3063
                      do i=1,molcrys%n_mol
3064
                         if (naa > molcrys%mol(i)%natoms) then
3065
                            naa=naa-molcrys%mol(i)%natoms
3066
                            cycle
3067
                         end if
3068

    
3069
                         if (molcrys%mol(i)%lbiso(naa) ==0) then
3070
                            np_refi=np_refi+1
3071
                            V_Vec(np_refi)=molcrys%mol(i)%biso(naa)
3072
                            V_Name(np_refi)=trim(code_nam(4))//trim(molcrys%mol(i)%AtName(naa))
3073
                            molcrys%mol(i)%mbiso(naa)=1.0
3074
                            molcrys%mol(i)%lbiso(naa)=np_refi
3075
                            V_Bounds(1,np_refi)=xl
3076
                            V_Bounds(2,np_refi)=xu
3077
                            V_Bounds(3,np_refi)=xs
3078
                            V_BCon(np_refi)=ic
3079
                            V_list(np_refi)=na
3080
                         end if
3081
                      end do
3082
                   end if
3083

    
3084
                   if (na <= molcrys%n_free) then
3085
                      do j=1,6
3086
                         if (molcrys%atm(na)%lu(j) ==0) then
3087
                            molcrys%atm(na)%mu(j)=1.0
3088
                            call get_atombet_ctr(molcrys%atm(na)%x,molcrys%atm(na)%u,molcrys%Spg, &
3089
                                                 np_refi,molcrys%atm(na)%lu,molcrys%atm(na)%mu)
3090
                            if (molcrys%atm(na)%lu(j) == np_refi) then
3091
                               V_Vec(np_refi)=molcrys%atm(na)%u(j)
3092
                               V_Name(np_refi)=trim(code_nam(5+j))//trim(molcrys%atm(na)%lab)
3093
                               V_Bounds(1,np_refi)=xl
3094
                               V_Bounds(2,np_refi)=xu
3095
                               V_Bounds(3,np_refi)=xs
3096
                               V_BCon(np_refi)=ic
3097
                               V_list(np_refi)=na
3098
                            else
3099
                               np_refi=np_refi-1
3100
                            end if
3101
                         end if
3102
                      end do
3103
                   end if
3104

    
3105
                   select case (nmol)
3106
                      case (-1)
3107
                         err_refcodes=.true.
3108
                         ERR_RefCodes_Mess="Option Not defined"
3109
                         return
3110

    
3111
                      case (0)
3112
                         do i=1,molcrys%n_mol
3113
                            write(unit=car,fmt="(i2)") i
3114
                            car=adjustl(car)
3115
                            do j=1,3
3116
                               if (molcrys%mol(i)%lxcentre(j) ==0) then
3117
                                  np_refi=np_refi+1
3118
                                  V_Vec(np_refi)=molcrys%mol(i)%xcentre(j)
3119
                                  V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3120
                                  molcrys%mol(i)%mxcentre(j)=1.0
3121
                                  molcrys%mol(i)%lxcentre(j)=np_refi
3122
                                  V_Bounds(1,np_refi)=xl
3123
                                  V_Bounds(2,np_refi)=xu
3124
                                  V_Bounds(3,np_refi)=xs
3125
                                  V_BCon(np_refi)=ic
3126
                                  V_list(np_refi)=-i
3127
                               end if
3128
                            end do
3129
                         end do
3130

    
3131
                      case (1:)
3132
                         write(unit=car,fmt="(i2)") nmol
3133
                         car=adjustl(car)
3134
                         do j=1,3
3135
                            if (molcrys%mol(nmol)%lxcentre(j) ==0) then
3136
                               np_refi=np_refi+1
3137
                               V_Vec(np_refi)=molcrys%mol(nmol)%xcentre(j)
3138
                               V_Name(np_refi)=trim(code_nam(12+j))//"entre_Mol"//trim(car)
3139
                               molcrys%mol(nmol)%mxcentre(j)=1.0
3140
                               molcrys%mol(nmol)%lxcentre(j)=np_refi
3141
                               V_Bounds(1,np_refi)=xl
3142
                               V_Bounds(2,np_refi)=xu
3143
                               V_Bounds(3,np_refi)=xs
3144
                               V_BCon(np_refi)=ic
3145
                               V_list(np_refi)=-nmol
3146
                            end if
3147
                         end do
3148
                   end select
3149

    
3150
                   select case (nmol)
3151
                      case (-1)
3152
                         err_refcodes=.true.
3153
                         ERR_RefCodes_Mess="Option Not defined"
3154
                         return
3155

    
3156
                      case (0)
3157
                         do i=1,molcrys%n_mol
3158
                            write(unit=car,fmt="(i2)") i
3159
                            car=adjustl(car)
3160
                            do j=1,3
3161
                               if (molcrys%mol(i)%lOrient(j) ==0) then
3162
                                  np_refi=np_refi+1
3163
                                  V_Vec(np_refi)=molcrys%mol(i)%Orient(j)
3164
                                  V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3165
                                  molcrys%mol(i)%mOrient(j)=1.0
3166
                                  molcrys%mol(i)%lOrient(j)=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)=-i
3172
                               end if
3173
                            end do
3174
                         end do
3175

    
3176
                      case (1:)
3177
                         write(unit=car,fmt="(i2)") nmol
3178
                         car=adjustl(car)
3179
                         do j=1,3
3180
                            if (molcrys%mol(nmol)%lOrient(j) ==0) then
3181
                               np_refi=np_refi+1
3182
                               V_Vec(np_refi)=molcrys%mol(nmol)%Orient(j)
3183
                               V_Name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"//trim(car)
3184
                               molcrys%mol(nmol)%mOrient(j)=1.0
3185
                               molcrys%mol(nmol)%lOrient(j)=np_refi
3186
                               V_Bounds(1,np_refi)=xl
3187
                               V_Bounds(2,np_refi)=xu
3188
                               V_Bounds(3,np_refi)=xs
3189
                               V_BCon(np_refi)=ic
3190
                               V_list(np_refi)=-nmol
3191
                            end if
3192
                         end do
3193
                   end select
3194

    
3195
                case (6)
3196
                   !---- CEN ----!
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
                case (7)
3243
                   !---- ORI  ----!
3244
                   select case (nmol)
3245
                      case (-1)
3246
                         err_refcodes=.true.
3247
                         ERR_RefCodes_Mess="Option Not defined"
3248
                         return
3249

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

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

    
3289
                case (8)
3290
                   !!! Not yet implemented !!!
3291

    
3292
             end select
3293
       end select
3294

    
3295
       return
3296
    End Subroutine Fill_RefCodes_Molcrys
3297

    
3298
    !!--++
3299
    !!--++ Subroutine Fill_RefCodes_Molec(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molec,Spg)
3300
    !!--++    integer,                      intent(in)     :: Key
3301
    !!--++    character(len=*),             intent(in)     :: Dire
3302
    !!--++    integer,                      intent(in)     :: Na
3303
    !!--++    integer,                      intent(in)     :: Nb
3304
    !!--++    real(kind=cp),                intent(in)     :: Xl
3305
    !!--++    real(kind=cp),                intent(in)     :: Xu
3306
    !!--++    real(kind=cp),                intent(in)     :: Xs
3307
    !!--++    integer,                      intent(in)     :: Ic
3308
    !!--++    type(molecule_type),          intent(in out) :: Molec
3309
    !!--++    type(space_group_type),       intent(in)     :: Spg
3310
    !!--++
3311
    !!--++ Overloaded
3312
    !!--++ Write on Vectors the Information for Molecule_Type
3313
    !!--++
3314
    !!--++ Update: March - 2005
3315
    !!
3316
    Subroutine Fill_RefCodes_Molec(Key,Dire,Na,Nb,Xl,Xu,Xs,Ic,Molec,Spg)
3317
       !---- Arguments ----!
3318
       integer,                      intent(in)     :: Key
3319
       character(len=*),             intent(in)     :: Dire
3320
       integer,                      intent(in)     :: Na
3321
       integer,                      intent(in)     :: Nb
3322
       real(kind=cp),                intent(in)     :: Xl
3323
       real(kind=cp),                intent(in)     :: Xu
3324
       real(kind=cp),                intent(in)     :: Xs
3325
       integer,                      intent(in)     :: Ic
3326
       type(molecule_type),          intent(in out) :: Molec
3327
       type(space_group_type),       intent(in)     :: Spg
3328

    
3329
       !---- Local variables ----!
3330
       integer :: j, nc
3331

    
3332
       call init_err_refcodes()
3333
       if (Na <= 0) then
3334
          err_refcodes=.true.
3335
          ERR_RefCodes_Mess="Number of atom no defined"
3336
          return
3337
       end if
3338

    
3339
       select case (dire)
3340
          !---- FIX Directive ----!
3341
          case ("fix")
3342

    
3343
             select case (key)
3344
                case (0)
3345

    
3346
                   !---- nb must be different zero ----!
3347
                   select case (nb)
3348
                      case (0)
3349
                         err_refcodes=.true.
3350
                         ERR_RefCodes_Mess="Option not defined"
3351
                         return
3352

    
3353
                      case ( 1:3)
3354
                         !---- X_, Y_, Z_----!
3355
                         if (molec%lI_coor(nb,na) /=0) then
3356
                            nc=molec%lI_coor(nb,na)
3357
                            call Delete_RefCodes(nc,molec)
3358
                         end if
3359

    
3360
                      case ( 4)
3361
                         !---- Biso_ ----!
3362
                         if (molec%lbiso(na) /=0) then
3363
                            nc=molec%lbiso(na)
3364
                            call Delete_RefCodes(nc,molec)
3365
                         end if
3366

    
3367
                      case ( 5)
3368
                         !---- Occ_ ----!
3369
                         if (molec%locc(na) /=0) then
3370
                            nc=molec%locc(na)
3371
                            call Delete_RefCodes(nc,molec)
3372
                         end if
3373

    
3374
                      case ( 6:12)
3375
                         !---- B11_, ..., B23_ ----!
3376
                         err_refcodes=.true.
3377
                         ERR_RefCodes_Mess="Option not defined"
3378
                         return
3379

    
3380
                      case (13:15)
3381
                         !---- Xc_, Yc_, Zc_ ----!
3382
                         if (molec%lxcentre(nb-12) /=0) then
3383
                            nc=molec%lxcentre(nb-12)
3384
                            call Delete_RefCodes(nc,molec)
3385
                         end if
3386

    
3387
                      case (16:18)
3388
                         !---- Theta_, Phi_, Chi_ ----!
3389
                         if (molec%lOrient(nb-15) /=0) then
3390
                            nc=molec%lOrient(nb-15)
3391
                            call Delete_RefCodes(nc,molec)
3392
                         end if
3393

    
3394
                      case (19:21)
3395
                         !!! Not yet implement !!!
3396

    
3397
                   end select ! nb
3398

    
3399
                case (1)
3400
                   !---- XYZ ----!
3401
                   do j=1,3
3402
                      if (molec%lI_coor(j,na) /=0) then
3403
                         nc=molec%lI_coor(j,na)
3404
                         call Delete_RefCodes(nc,molec)
3405
                      end if
3406
                   end do
3407

    
3408
                case (2)
3409
                   !---- OCC ----!
3410
                   if (molec%locc(na) /=0) then
3411
                      nc=molec%locc(na)
3412
                      call Delete_RefCodes(nc,molec)
3413
                   end if
3414

    
3415
                case (3)
3416
                   !---- BIS ----!
3417
                   if (molec%lbiso(na) /=0) then
3418
                      nc=molec%lbiso(na)
3419
                      call Delete_RefCodes(nc,molec)
3420
                   end if
3421

    
3422
                case (4)
3423
                   !---- BAN ----!
3424
                   err_refcodes=.true.
3425
                   ERR_RefCodes_Mess="Option not defined"
3426
                   return
3427

    
3428
                case (5)
3429
                  !---- ALL ----!
3430
                  do j=1,3
3431
                     if (molec%lI_coor(j,na) /=0) then
3432
                        nc=molec%lI_coor(j,na)
3433
                        call Delete_RefCodes(nc,molec)
3434
                     end if
3435
                  end do
3436
                  if (molec%locc(na) /=0) then
3437
                     nc=molec%locc(na)
3438
                     call Delete_RefCodes(nc,molec)
3439
                  end if
3440
                  if (molec%lbiso(na) /=0) then
3441
                     nc=molec%lbiso(na)
3442
                     call Delete_RefCodes(nc,molec)
3443
                  end if
3444
                  do j=1,3
3445
                     if (molec%lxcentre(j) /=0) then
3446
                        nc=molec%lxcentre(j)
3447
                        call Delete_RefCodes(nc,molec)
3448
                     end if
3449
                  end do
3450
                  do j=1,3
3451
                     if (molec%lOrient(j) /=0) then
3452
                        nc=molec%lOrient(j)
3453
                        call Delete_RefCodes(nc,molec)
3454
                     end if
3455
                  end do
3456
                  !!! Falta THermal TLS !!!
3457

    
3458
               case (6)
3459
                  !---- CEN ----!
3460
                  do j=1,3
3461
                     if (molec%lxcentre(j) /=0) then
3462
                        nc=molec%lxcentre(j)
3463
                        call Delete_RefCodes(nc,molec)
3464
                     end if
3465
                  end do
3466

    
3467
               case (7)
3468
                  !---- ORI ----!
3469
                  do j=1,3
3470
                     if (molec%lOrient(j) /=0) then
3471
                        nc=molec%lOrient(j)
3472
                        call Delete_RefCodes(nc,molec)
3473
                     end if
3474
                  end do
3475

    
3476
               case (8)
3477
                  !---- THE ----!
3478
                  !!! Not Yet Implemented !!!
3479
             end select
3480

    
3481
          !---- VARY Directive ----!
3482
          case ("var")
3483

    
3484
             select case (key)
3485
                case (0)
3486

    
3487
                   !---- Nb must be different zero ----!
3488
                   select case (nb)
3489
                      case (0)
3490
                         err_refcodes=.true.
3491
                         ERR_RefCodes_Mess="Option not defined"
3492
                         return
3493

    
3494
                      case ( 1:3)
3495
                         !--- X_, Y_, Z_ ----!
3496
                         if (molec%lI_coor(nb,na) ==0) then
3497
                            molec%mI_coor(nb,na)=1.0
3498
                            call get_atompos_ctr(molec%I_Coor(:,na),  &
3499
                                                 Spg, np_refi,  &
3500
                                                 molec%lI_coor(:,na), &
3501
                                                 molec%mI_coor(:,na))
3502
                            if (molec%lI_coor(nb,na) == np_refi) then
3503
                               V_Vec(np_refi)=molec%I_Coor(nb,na)
3504
                               V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3505
                               V_Bounds(1,np_refi)=xl
3506
                               V_Bounds(2,np_refi)=xu
3507
                               V_Bounds(3,np_refi)=xs
3508
                               V_BCon(np_refi)=ic
3509
                               V_list(np_refi)=na
3510
                            else
3511
                               np_refi=np_refi-1
3512
                            end if
3513
                         end if
3514

    
3515
                      case ( 4)
3516
                         !---- Biso_ ----!
3517
                         if (molec%lbiso(na) ==0) then
3518
                            np_refi=np_refi+1
3519
                            V_Vec(np_refi)=molec%biso(na)
3520
                            V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3521
                            molec%mbiso(na)=1.0
3522
                            molec%lbiso(na)=np_refi
3523
                            V_Bounds(1,np_refi)=xl
3524
                            V_Bounds(2,np_refi)=xu
3525
                            V_Bounds(3,np_refi)=xs
3526
                            V_BCon(np_refi)=ic
3527
                            V_list(np_refi)=na
3528
                         end if
3529

    
3530
                      case ( 5)
3531
                         !---- Occ_ ----!
3532
                         if (molec%locc(na) ==0) then
3533
                            np_refi=np_refi+1
3534
                            V_Vec(np_refi)=molec%occ(na)
3535
                            V_Name(np_refi)=trim(code_nam(nb))//trim(molec%AtName(na))
3536
                            molec%mocc(na)=1.0
3537
                            molec%locc(na)=np_refi
3538
                            V_Bounds(1,np_refi)=xl
3539
                            V_Bounds(2,np_refi)=xu
3540
                            V_Bounds(3,np_refi)=xs
3541
                            V_BCon(np_refi)=ic
3542
                            V_list(np_refi)=na
3543
                         end if
3544

    
3545
                      case ( 6:12)
3546
                         !---- B11_, ..., B23_ ----!
3547
                         err_refcodes=.true.
3548
                         ERR_RefCodes_Mess="Option not defined"
3549
                         return
3550

    
3551
                      case (13:15)
3552
                         !---- Xc_, Yc_, Zc_ ----!
3553
                         if (molec%lxcentre(nb-12) ==0) then
3554
                            np_refi=np_refi+1
3555
                            V_Vec(np_refi)=molec%xcentre(nb-12)
3556
                            V_Name(np_refi)=trim(code_nam(nb))//"Mol"
3557
                            molec%mxcentre(nb-12)=1.0
3558
                            molec%lxcentre(nb-12)=np_refi
3559
                            V_Bounds(1,np_refi)=xl
3560
                            V_Bounds(2,np_refi)=xu
3561
                            V_Bounds(3,np_refi)=xs
3562
                            V_BCon(np_refi)=ic
3563
                            V_list(np_refi)=0
3564
                         end if
3565

    
3566
                      case (16:18)
3567
                         !---- Theta_, Phi_, Chi_ ----!
3568
                         if (molec%lOrient(nb-15) ==0) then
3569
                            np_refi=np_refi+1
3570
                            V_Vec(np_refi)=molec%orient(nb-15)
3571
                            V_Name(np_refi)=trim(code_nam(nb))//"Mol"
3572
                            molec%mOrient(nb-15)=1.0
3573
                            molec%lOrient(nb-15)=np_refi
3574
                            V_Bounds(1,np_refi)=xl
3575
                            V_Bounds(2,np_refi)=xu
3576
                            V_Bounds(3,np_refi)=xs
3577
                            V_BCon(np_refi)=0
3578
                            V_list(np_refi)=0
3579
                         end if
3580

    
3581
                      case (19:21)
3582
                         !!! Not yet implement !!!
3583

    
3584
                   end select ! nb
3585

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

    
3609
                case (2)
3610
                   !---- OCC ----!
3611
                   if (molec%locc(na) ==0) then
3612
                      np_refi=np_refi+1
3613
                      V_Vec(np_refi)=molec%occ(na)
3614
                      V_Name(np_refi)=trim(code_nam(5))//trim(molec%AtName(na))
3615
                      molec%mocc(na)=1.0
3616
                      molec%locc(na)=np_refi
3617
                      V_Bounds(1,np_refi)=xl
3618
                      V_Bounds(2,np_refi)=xu
3619
                      V_Bounds(3,np_refi)=xs
3620
                      V_BCon(np_refi)=ic
3621
                      V_list(np_refi)=na
3622
                   end if
3623

    
3624
                case (3)
3625
                   !---- BIS ----!
3626
                   if (molec%lbiso(na) ==0) then
3627
                      np_refi=np_refi+1
3628
                      V_Vec(np_refi)=molec%biso(na)
3629
                      V_Name(np_refi)=trim(code_nam(4))//trim(molec%AtName(na))
3630
                      molec%mbiso(na)=1.0
3631
                      molec%lbiso(na)=np_refi
3632
                      V_Bounds(1,np_refi)=xl
3633
                      V_Bounds(2,np_refi)=xu
3634
                      V_Bounds(3,np_refi)=xs
3635
                      V_BCon(np_refi)=ic
3636
                      V_list(np_refi)=na
3637
                   end if
3638

    
3639
                case (4)
3640
                   !---- BAN ----!
3641
                   err_refcodes=.true.
3642
                   ERR_RefCodes_Mess="Option not defined"
3643
                   return
3644

    
3645
                case (5)
3646
                   !---- ALL ----!
3647
                   do j=1,3
3648
                      if (molec%lI_coor(j,na) ==0) then
3649
                         molec%mI_coor(j,na)=1.0
3650
                         call get_atompos_ctr(molec%I_Coor(:,na),  &
3651
                                              Spg, np_refi,  &
3652
                                              molec%lI_coor(:,na), &
3653
                                              molec%mI_coor(:,na))
3654
                         if (molec%lI_coor(j,na) == np_refi) then
3655
                            V_Vec(np_refi)=molec%I_Coor(j,na)
3656
                            V_Name(np_refi)=trim(code_nam(j))//trim(molec%AtName(na))
3657
                            V_Bounds(1,np_refi)=xl
3658
                            V_Bounds(2,np_refi)=xu
3659
                            V_Bounds(3,np_refi)=xs
3660
                            V_BCon(np_refi)=ic
3661
                            V_list(np_refi)=na
3662
                         else
3663
                            np_refi=np_refi-1
3664
                         end if
3665
                      end if
3666
                   end do
3667
                   if (molec%locc(na) ==0) then
3668
                      np_refi=np_refi+1
3669
                      V_Vec(np_refi)=molec%occ(na)
3670
                      V_Name(np_refi)=trim(code_nam(5))//trim(molec%AtName(na))
3671
                      molec%mocc(na)=1.0
3672
                      molec%locc(na)=np_refi
3673
                      V_Bounds(1,np_refi)=xl
3674
                      V_Bounds(2,np_refi)=xu
3675
                      V_Bounds(3,np_refi)=xs
3676
                      V_BCon(np_refi)=ic
3677
                      V_list(np_refi)=na
3678
                   end if
3679
                   if (molec%lbiso(na) ==0) then
3680
                      np_refi=np_refi+1
3681
                      V_Vec(np_refi)=molec%biso(na)
3682
                      V_Name(np_refi)=trim(code_nam(4))//trim(molec%AtName(na))
3683
                      molec%mbiso(na)=1.0
3684
                      molec%lbiso(na)=np_refi
3685
                      V_Bounds(1,np_refi)=xl
3686
                      V_Bounds(2,np_refi)=xu
3687
                      V_Bounds(3,np_refi)=xs
3688
                      V_BCon(np_refi)=ic
3689
                      V_list(np_refi)=na
3690
                   end if
3691
                   do j=1,3
3692
                      if (molec%lxcentre(j) ==0) then
3693
                         np_refi=np_refi+1
3694
                         V_Vec(np_refi)=molec%xcentre(j)
3695
                         V_name(np_refi)=trim(code_nam(12+j))//"entre_Mol"
3696
                         molec%mxcentre(j)=1.0
3697
                         molec%lxcentre(j)=np_refi
3698
                         V_Bounds(1,np_refi)=xl
3699
                         V_Bounds(2,np_refi)=xu
3700
                         V_Bounds(3,np_refi)=xs
3701
                         V_BCon(np_refi)=ic
3702
                         V_list(np_refi)=0
3703
                      end if
3704
                   end do
3705
                   do j=1,3
3706
                      if (molec%lOrient(j) ==0) then
3707
                         np_refi=np_refi+1
3708
                         V_Vec(np_refi)=molec%orient(j)
3709
                         V_name(np_refi)=trim(code_nam(15+j))//"Orient_Mol"
3710
                         molec%mOrient(j)=1.0
3711
                         molec%lOrient(j)=np_refi
3712
                         V_Bounds(1,np_refi)=xl
3713
                         V_Bounds(2,np_refi)=xu
3714
                         V_Bounds(3,np_refi)=xs
3715
                         V_BCon(np_refi)=0
3716
                         V_list(np_refi)=0
3717
                      end if
3718
                   end do
3719

    
3720
                   !!! Falta THE !!!
3721

    
3722
                case (6)
3723
                   !---- CEN ----!
3724
                   do j=1,3
3725
                      if (molec%lxcentre(j) ==0) then
3726
                         np_refi=np_refi+1
3727
                         V_Vec(np_refi)=molec%xcentre(j)
3728
                         V_name(np_refi)=trim(code_nam(12+j))//"Mol"
3729
                         molec%mxcentre(j)=1.0
3730
                         molec%lxcentre(j)=np_refi
3731
                         V_Bounds(1,np_refi)=xl
3732
                         V_Bounds(2,np_refi)=xu
3733
                         V_Bounds(3,np_refi)=xs
3734
                         V_BCon(np_refi)=ic
3735
                         V_list(np_refi)=0
3736
                      end if
3737
                   end do
3738

    
3739
                case (7)
3740
                   !---- ORI ----!
3741
                   do j=1,3
3742
                      if (molec%lOrient(j) ==0) then
3743
                         np_refi=np_refi+1
3744
                         V_Vec(np_refi)=molec%orient(j)
3745
                         V_name(np_refi)=trim(code_nam(15+j))//"Mol"
3746
                         molec%mOrient(j)=1.0
3747
                         molec%lOrient(j)=np_refi
3748
                         V_Bounds(1,np_refi)=xl
3749
                         V_Bounds(2,np_refi)=xu
3750
                         V_Bounds(3,np_refi)=xs
3751
                         V_BCon(np_refi)=ic
3752
                         V_list(np_refi)=0
3753
                      end if
3754
                   end do
3755

    
3756
                case (8)
3757
                   !---- THE ----!
3758

    
3759
                   !!! Not yet implemented !!!
3760
             end select
3761
       end select
3762

    
3763
       return
3764
    End Subroutine Fill_RefCodes_Molec
3765

    
3766
    !!--++
3767
    !!--++  Subroutine Get_Atombet_Ctr(X,Betas,Spgr,Codini,Icodes,Multip,Ord,Ss,Ipr)
3768
    !!--++     real(kind=cp), dimension(3),             intent(in    ) :: X         !Atom position (fractional coordinates)
3769
    !!--++     real(kind=cp), dimension(6),             intent(in out) :: Betas     !Anisotropic temperature factors
3770
    !!--++     type(Space_Group_type),                  intent(in    ) :: Spgr      !Space Group
3771
    !!--++     Integer,                                 intent(in out) :: Codini    !Last attributed parameter
3772
    !!--++     Integer, dimension(6),                   intent(in out) :: Icodes    !codewords for betas only number
3773
    !!--++     real(kind=cp), dimension(6),             intent(in out) :: Multip    !Multipliers
3774
    !!--++     integer,                       optional, intent(in    ) :: Ord       !Order of the stabilizer
3775
    !!--++     integer, dimension(:),         optional, intent(in    ) :: Ss        !Pointer to SymmOp. of stabilizer
3776
    !!--++     integer,                       optional, intent(in    ) :: Ipr       !Printing unit for debug
3777
    !!--++
3778
    !!--++  Subroutine to get the appropriate constraints in the refinement codes of
3779
    !!--++  anisotropic atomic displacement(thermal) parameters.
3780
    !!--++  New algorithm based in the Wigner theorem.
3781
    !!--++  The matrix Bet = Sum { R Beta RT} displays the symmetry constraints to be
3782
    !!--++  applied to the anisotropic temperature factors. The sum runs over all rotational
3783
    !!--++  symmetry operators of the stabilizer of the particular atom position in the given
3784
    !!--++  space group.
3785
    !!--++  There are a total of 29 kind of relations that may appear in the Bet matrix:
3786
    !!--++
3787
    !!--++     1    A A A 0   0   0  -> m-3m, -43m, 432, m-3,23, 3[111].2[001]
3788
    !!--++     2    A A C 0   0   0  -> 4/mmm, -42m, 4mm, 422, 4/m, -4,4, 4[001]
3789
    !!--++     3    A B A 0   0   0  -> 4[010]
3790
    !!--++     4    A B B 0   0   0  -> 4[100]
3791
    !!--++     5    A A A D   D   D  -> -3m, 3m, 32, -3, 3   3[111]
3792
    !!--++     6    A A A D  -D  -D  -> 3[11-1]
3793
    !!--++     7    A A A D  -D   D  -> 3[1-11]
3794
    !!--++     8    A A A D   D  -D  -> 3[-111]
3795
    !!--++     9    A A C A/2 0   0  -> 6/mmm, -6m2, 6mm, 622, 6/m, 6,-6,-3m, 32,-3, 3:  h 3[001]
3796
    !!--++    10    A B C 0   0   0  -> mmm, mm2, 222  2[001] 2[100]
3797
    !!--++    11    A A C D   0   0  -> 2[001], 2[110]    w
3798
    !!--++    12    A B A 0   E   0  -> 2[010], 2[101]
3799
    !!--++    13    A B B 0   0   F  -> 2[100], 2[011]
3800
    !!--++    14    A B C B/2 0   0  -> 2[001], 2[100]    h
3801
    !!--++    15    A B C A/2 0   0  -> 2[001], 2[010]    h
3802
    !!--++    16    A B C D   0   0  -> 2/m, m, 2: 2[001] w
3803
    !!--++    17    A B C 0   E   0  -> 2[010]
3804
    !!--++    18    A B C 0   0   F  -> 2[100]
3805
    !!--++    19    A A C D   E  -E  -> 2[110]            w
3806
    !!--++    20    A A C D   E   E  -> 2[1-10]           w
3807
    !!--++    21    A B A D   E  -D  -> 2[101]
3808
    !!--++    22    A B A D   E   D  -> 2[10-1]
3809
    !!--++    23    A B B D  -D   F  -> 2[011]
3810
    !!--++    24    A B B D   D   F  -> 2[01-1]
3811
    !!--++    25    A B C B/2 F/2 F  -> 2[100]            h
3812
    !!--++    26    A B C A/2 0   F  -> 2[210]            h
3813
    !!--++    27    A B C B/2 E   0  -> 2[120]            h
3814
    !!--++    28    A B C A/2 E   E/2-> 2[010]            h
3815
    !!--++    29    A B C D   E   F  -> 1, -1
3816
    !!--++
3817
    !!--++   Updated: 14 July 2011
3818
    !!--++
3819
    !!
3820

    
3821
    Subroutine Get_Atombet_Ctr(X,Betas,Spgr,Codini,Icodes,Multip,Ord,Ss,Ipr)
3822
       real(kind=cp), dimension(3),             intent(in    ) :: X
3823
       real(kind=cp), dimension(6),             intent(in out) :: Betas
3824
       type(Space_Group_type),                  intent(in    ) :: Spgr
3825
       Integer,                                 intent(in out) :: Codini
3826
       Integer, dimension(6),                   intent(in out) :: Icodes
3827
       real(kind=cp), dimension(6),             intent(in out) :: Multip
3828
       integer,                       optional, intent(in    ) :: Ord
3829
       integer, dimension(:),         optional, intent(in    ) :: Ss
3830
       integer,                       optional, intent(in    ) :: Ipr
3831

    
3832

    
3833
       ! Local variables
3834
       character (len=1), dimension(6)   :: cdd
3835
       integer                           :: i,j,order
3836
       integer,           dimension(48)  :: ss_ptr
3837
       integer,           dimension(6)   :: codd
3838
       integer,           dimension(3,3) :: Rsym
3839
       real(kind=cp),     dimension(3,3) :: bet,bett,Rs
3840
       real(kind=cp),     dimension(6)   :: cod
3841
       real(kind=cp),     dimension(3,48):: atr
3842
       real(kind=cp),     parameter      :: epss=0.01_cp
3843

    
3844
       cod=real(icodes)
3845

    
3846
       do j=1,6
3847
          if (cod(j) < 1.0 .and. abs(multip(j)) > epss)  then
3848
             codini=codini+1
3849
             cod(j) = real(codini)
3850
          end if
3851
       end do
3852

    
3853
       if (present(ord) .and. present(ss)) then
3854
          order=ord
3855
          ss_ptr(1:order) = ss(1:ord)
3856
       else
3857
          call get_stabilizer(x,Spgr,order,ss_ptr,atr)
3858
       end if
3859

    
3860
       bet=reshape((/17.0, 7.0,3.0,  &
3861
                     7.0,13.0,5.0,  &
3862
                     3.0, 5.0,11.0/),(/3,3/))
3863
       bett=bet
3864
       if (order > 1 ) then
3865
          do j=2,order
3866
             Rsym=Spgr%SymOp(ss_ptr(j))%Rot
3867
             Rs=real(Rsym)
3868
             bett=bett+ matmul(Rs,matmul(bet,transpose(Rs)))
3869
          end do
3870
       end if
3871
       Rsym=nint(1000.0*bett)
3872
       codd=(/Rsym(1,1),Rsym(2,2),Rsym(3,3),Rsym(1,2),Rsym(1,3),Rsym(2,3)/)
3873
       cdd=(/'a','b','c','d','e','f'/)
3874
       multip=1.0
3875
       !Search systematically all the possible constraints
3876

    
3877
       if(codd(1) == codd(2) .and. codd(1) == codd(3)) then ! a a a
3878
         if(codd(4) == codd(5) .and. codd(4) == codd(6) ) then ! a a a d d d
3879
           if(codd(4) == 0) then
3880
             cdd=(/'a','a','a','0','0','0'/)     ! 1 A A A 0   0   0
3881
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
3882
             betas(4:6)=0.0
3883
             betas(2:3)=betas(1)
3884
             cod(2:3)=cod(1); cod(4:6)=0.0
3885
           else
3886
             cdd=(/'a','a','a','d','d','d'/)     ! 5 A A A D   D   D
3887
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
3888
             betas(5:6)=betas(4)
3889
             betas(2:3)=betas(1)
3890
             cod(2:3)=cod(1); cod(5:6)=cod(4)
3891
           end if
3892
         else if(codd(4) == -codd(5) .and. codd(4) == -codd(6) ) then !a a a d -d -d
3893
           cdd=(/'a','a','a','d','d','d'/)       ! 6 A A A D  -D  -D
3894
           multip=(/1.0,1.0,1.0,1.0,-1.0,-1.0/)
3895
           betas(5:6)=-betas(4)
3896
           betas(2:3)=betas(1)
3897
           cod(2:3)=cod(1); cod(5:6)=cod(4)
3898
         else if(codd(4) == -codd(5) .and. codd(4) ==  codd(6) ) then !a a a d -d  d
3899
           cdd=(/'a','a','a','d','d','d'/)       ! 7 A A A D  -D   D
3900
           multip=(/1.0,1.0,1.0,1.0,-1.0, 1.0/)
3901
           betas(5)=-betas(4); betas(6)=betas(4)
3902
           betas(2:3)=betas(1)
3903
           cod(2:3)=cod(1); cod(5:6)= cod(4)
3904
         else if(codd(4) ==  codd(5) .and. codd(4) == -codd(6) ) then !a a a d  d -d
3905
           cdd=(/'a','a','a','d','d','d'/)       ! 8 A A A D   D  -D
3906
           multip=(/1.0,1.0,1.0,1.0, 1.0,-1.0/)
3907
           betas(6)=-betas(4); betas(5)=betas(4)
3908
           betas(2:3)=betas(1)
3909
           cod(2:3)=cod(1); cod(5:6)= cod(4)
3910
         end if
3911

    
3912
       else if(codd(1) == codd(2)) then ! a a c
3913
         if(codd(4) == codd(5) .and. codd(4) == codd(6) .and. codd(4) == 0) then ! a a c 0 0 0
3914
             cdd=(/'a','a','c','0','0','0'/)     ! 2 A A C 0   0   0
3915
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
3916
             betas(4:6)=0.0
3917
             betas(2)=betas(1)
3918
             cod(2)=cod(1); cod(4:6)= 0.0
3919
         else if(codd(5) == codd(6) .and. codd(5) == 0) then ! a a c x 0 0
3920
             if(codd(4) == codd(1)/2) then
3921
               cdd=(/'a','a','c','a','0','0'/)     ! 9 A A C A/2 0   0
3922
               multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
3923
               betas(5:6)=0.0; betas(4)=betas(1)*0.5
3924
               betas(2)=betas(1)
3925
               cod(2)=cod(1); cod(4)= cod(1); cod(5:6)=0.0
3926
             else
3927
               cdd=(/'a','a','c','d','0','0'/)     !11 A A C D   0   0
3928
               multip=(/1.0,1.0,1.0,1.0,0.0,0.0/)
3929
               betas(5:6)=0.0
3930
               betas(2)=betas(1)
3931
               cod(2)=cod(1); cod(5:6)=0.0
3932
             end if
3933
         else
3934
             if(codd(5) == codd(6)) then  ! a a c d e e
3935
               cdd=(/'a','a','c','d','e','e'/)     !20 A A C D   E   E
3936
               multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
3937
               betas(6)=betas(5)
3938
               betas(2)=betas(1)
3939
               cod(2)=cod(1); cod(6)=cod(5)
3940
             else if(codd(5) == -codd(6)) then  ! a a c d e -e
3941
               cdd=(/'a','a','c','d','e','e'/)     !19 A A C D   E  -E
3942
               multip=(/1.0,1.0,1.0,1.0,1.0,-1.0/)
3943
               betas(6)=-betas(5)
3944
               betas(2)=betas(1)
3945
               cod(2)=cod(1); cod(6)=cod(5)
3946
             end if
3947
         end if
3948

    
3949
       else if(codd(1) == codd(3)) then ! a b a
3950
         if(codd(4) == codd(6)) then    ! a b a d x d
3951
           if(codd(4) == 0) then  ! a b a 0 x 0
3952
             if(codd(5) == 0) then ! a b a 0 0 0
3953
               cdd=(/'a','b','a','0','0','0'/)     ! 3 A B A 0   0   0
3954
               multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
3955
               betas(4:6)=0.0
3956
               betas(3)=betas(1)
3957
               cod(3)=cod(1); cod(4:6)=0.0
3958
             else                  ! a b a 0 e 0
3959
               cdd=(/'a','b','a','0','e','0'/)     !12 A B A 0   E   0
3960
               multip=(/1.0,1.0,1.0,0.0,1.0,0.0/)
3961
               betas(4)=0.0;  betas(6)=0.0
3962
               betas(3)=betas(1)
3963
               cod(3)=cod(1); cod(4)=0.0;  cod(6)=0.0
3964
             endif
3965
           else  !! a b a d e d
3966
             cdd=(/'a','b','a','d','e','d'/)       !22 A B A D   E   D
3967
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
3968
             betas(6)=betas(4)
3969
             betas(3)=betas(1)
3970
             cod(3)=cod(1); cod(6)=cod(4)
3971
          end if
3972

    
3973
         else if(codd(4) == -codd(6)) then ! a b a d e -d
3974
           cdd=(/'a','b','a','d','e','d'/)         !21 A B A D   E  -D
3975
           multip=(/1.0,1.0,1.0,1.0,1.0,-1.0/)
3976
           betas(6)=-betas(4)
3977
           betas(3)=betas(1)
3978
           cod(3)=cod(1); cod(6)=cod(4)
3979
         end if
3980

    
3981
       else if(codd(2) == codd(3)) then ! a b b
3982
         if(codd(4) == codd(5)) then    ! a b b d d x
3983
           if(codd(4) == 0) then  ! a b b 0 0 x
3984
             if(codd(6) == 0) then ! a b b 0 0 0
3985
               cdd=(/'a','b','b','0','0','0'/)     ! 4 A B B 0   0   0
3986
               multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
3987
               betas(4:6)=0.0
3988
               betas(3)=betas(2)
3989
               cod(3)=cod(2); cod(4:6)=0.0
3990
             else                  ! a b b 0 0 f
3991
               cdd=(/'a','b','b','0','0','f'/)     !13 A B B 0   0   F
3992
               multip=(/1.0,1.0,1.0,0.0,0.0,1.0/)
3993
               betas(4:5)=0.0
3994
               betas(3)=betas(2)
3995
               cod(3)=cod(2); cod(4:5)=0.0
3996
             endif
3997
           else  !! a b b d d f
3998
             cdd=(/'a','b','b','d','d','f'/)       !24 A B B D   D   F
3999
             multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4000
             betas(5)=betas(4)
4001
             betas(3)=betas(2)
4002
             cod(3)=cod(2); cod(5)=cod(4)
4003
           end if
4004
         else if(codd(4) == -codd(5)) then ! a b b d -d e
4005
           cdd=(/'a','b','b','d','d','f'/)         !23 A B B D  -D   F
4006
           multip=(/1.0,1.0,1.0,1.0,-1.0,1.0/)
4007
           betas(5)=-betas(4)
4008
           betas(3)=betas(2)
4009
           cod(3)=cod(2); cod(5)=cod(4)
4010
         end if
4011

    
4012
       else !Now a /= b /= c
4013

    
4014
         if(codd(4) == codd(5) .and. codd(4) == 0) then ! a b c 0 0 x
4015
           if(codd(6) == 0) then ! a b c 0 0 0
4016
             cdd=(/'a','b','c','0','0','0'/)          !10 A B C 0   0   0
4017
             multip=(/1.0,1.0,1.0,0.0,0.0,0.0/)
4018
             betas(4:6)=0.0
4019
             cod(4:6)=0.0
4020
           else
4021
             cdd=(/'a','b','c','0','0','f'/)          !18 A B C 0   0   F
4022
             multip=(/1.0,1.0,1.0,0.0,0.0,1.0/)
4023
             betas(4:5)=0.0
4024
             cod(4:5)=0.0
4025
           end  if
4026
         else if(codd(5) == codd(6) .and. codd(5) == 0) then  ! a b c x 0 0
4027
           if(codd(4) == codd(1)/2) then ! a b c a/2 0 0
4028
             cdd=(/'a','b','c','a','0','0'/)          !15 A B C A/2 0   0
4029
             multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
4030
             betas(5:6)=0.0; betas(4)=betas(1)*0.5
4031
             cod(4)=cod(1); cod(5:6)=0.0
4032
           else if(codd(4) == codd(2)/2) then    !a b c b/2 0 0
4033
             cdd=(/'a','b','c','b','0','0'/)          !14 A B C B/2 0   0
4034
             multip=(/1.0,1.0,1.0,0.5,0.0,0.0/)
4035
             betas(5:6)=0.0; betas(4)=betas(2)*0.5
4036
             cod(4)=cod(2); cod(5:6)=0.0
4037
           else
4038
             cdd=(/'a','b','c','d','0','0'/)          !16 A B C D   0   0
4039
             multip=(/1.0,1.0,1.0,1.0,0.0,0.0/)
4040
             betas(5:6)=0.0
4041
             cod(5:6)=0.0
4042
           end  if
4043
         else if(codd(4) == codd(6) .and. codd(4) == 0) then !a b c 0 e 0
4044
           cdd=(/'a','b','c','0','e','0'/)            !17 A B C 0   E   0
4045
           multip=(/1.0,1.0,1.0,0.0,1.0,0.0/)
4046
           betas(4)=0.0; betas(6)=0.0
4047
           cod(4)=0.0; cod(6)=0.0
4048
         else if(codd(4) == codd(1)/2 .and. codd(5) == 0) then !a b c a/2 0 f
4049
           cdd=(/'a','b','c','a','0','f'/)            !26 A B C A/2 0   F
4050
           multip=(/1.0,1.0,1.0,0.5,0.0,1.0/)
4051
           betas(4)=betas(1)*0.5; betas(5)=0.0
4052
           cod(4)=cod(1); cod(5)=0.0
4053
         else if(codd(4) == codd(2)/2 .and. codd(6) == 0) then !a b c b/2 e 0
4054
           cdd=(/'a','b','c','b','e','0'/)            !27 A B C B/2 E   0
4055
           multip=(/1.0,1.0,1.0,0.5,1.0,0.0/)
4056
           betas(4)=betas(2)*0.5; betas(6)=0.0
4057
           cod(4)=cod(2); cod(6)=0.0
4058
         else if(codd(4) == codd(2)/2 .and. codd(5) == codd(6)/2) then !a b c b/2 f/2 f
4059
           cdd=(/'a','b','c','b','f','f'/)            !25 A B C B/2 F/2 F
4060
           multip=(/1.0,1.0,1.0,0.5,0.5,1.0/)
4061
           betas(4)=betas(2)*0.5; betas(5)=betas(6)*0.5
4062
           cod(4)=cod(2); cod(5)=cod(6)
4063
         else if(codd(4) == codd(1)/2 .and. codd(6) == codd(5)/2) then !a b c a/2 e e/2
4064
           cdd=(/'a','b','c','a','e','e'/)            !28 A B C A/2 E   E/2
4065
           multip=(/1.0,1.0,1.0,0.5,1.0,0.5/)
4066
           betas(4)=betas(1)*0.5; betas(6)=betas(5)*0.5
4067
           cod(4)=cod(1); cod(6)=cod(5)
4068
         else
4069
           cdd=(/'a','b','c','d','e','f'/)            !29 A B C D   E   F
4070
           multip=(/1.0,1.0,1.0,1.0,1.0,1.0/)
4071
         end if
4072
       end if
4073

    
4074
       do j=1,6
4075
          if (multip(j) < epss .or. cdd(j) == "0" ) then
4076
             icodes(j) = 0
4077
          else
4078
             icodes(j) = nint(cod(j))
4079
          end if
4080
       end do
4081

    
4082
       if(present(Ipr)) then
4083
         Write(Ipr,'(a,6i5)')           '     Codes on Betas       :  ',Icodes
4084
         Write(Ipr,'(a,6(a,1x),6f7.3)') '     Codes and multipliers:  ',cdd,multip
4085
         Write(Ipr,'(a)')               '     Beta_TOT matrix:  '
4086
         Do I=1,3
4087
          Write(Ipr,'(a,3f12.4)')       '                      ',bett(i,:)
4088
         End Do
4089
       end if
4090
       return
4091
    End Subroutine Get_Atombet_Ctr
4092

    
4093
    !!--++
4094
    !!--++  Subroutine Get_Atompos_Ctr(X,Spgr,Codini,Icodes,Multip,Ord,Ss,Att,Ipr)
4095
    !!--++     real(kind=cp), dimension(3),       intent(in    ) :: x      !Atom position (fractional coordinates)
4096
    !!--++     type(Space_Group_type),            intent(in    ) :: Spgr   !Space Group
4097
    !!--++     Integer,                           intent(in out) :: codini !Last attributed parameter
4098
    !!--++     integer,       dimension(3),       intent(in out) :: Icodes
4099
    !!--++     real(kind=cp), dimension(3),       intent(in out) :: Multip
4100
    !!--++     integer,                 optional, intent(in    ) :: Ord
4101
    !!--++     integer, dimension(:),   optional, intent(in    ) :: Ss
4102
    !!--++     integer, dimension(:,:), optional, intent(in    ) :: Atr
4103
    !!--++     integer,                 optional, intent(in    ) :: Ipr
4104
    !!--++
4105
    !!--++     (Private)
4106
    !!--++     Subroutine to get the appropriate constraints in the refinement codes of
4107
    !!--++     atoms positions. The algorithm is based in an analysis of the symbol generated
4108
    !!--++     for the symmetry elements of the operators belonging to the stabilizer of the
4109
    !!--++     atom position. This routine operates with splitted codes in the sense of
4110
    !!--++     FullProf rules
4111
    !!--++
4112
    !!--++     Updated: July - 2011 (JRC, the old subroutine has been completely changed)
4113
    !!
4114
    Subroutine Get_Atompos_Ctr(X,Spgr,Codini,ICodes,Multip,Ord,Ss,Att,Ipr)
4115
       real(kind=cp), dimension(3),            intent(in)     :: X
4116
       type(Space_Group_type),                 intent(in)     :: Spgr
4117
       Integer,                                intent(in out) :: Codini
4118
       Integer,       dimension(3),            intent(in out) :: ICodes
4119
       real(kind=cp), dimension(3),            intent(in out) :: Multip
4120
       integer,                       optional,intent(in)     :: Ord
4121
       integer, dimension(:),         optional,intent(in)     :: Ss
4122
       real(kind=cp), dimension(:,:), optional,intent(in)     :: Att
4123
       integer,                       optional,intent(in)     :: Ipr
4124

    
4125
       ! Local variables
4126
       integer                          :: i,j,k,order,L,L1,L2,ipar,j1
4127
       integer,          dimension(3,3) :: RSym
4128
       integer,          dimension(48)  :: ss_ptr
4129
       real(kind=cp),    dimension(3,48):: atr
4130
       real(kind=cp),    dimension(3)   :: tr
4131

    
4132
       character(len=40)                :: symbol,tsymbol,sym_symb
4133
       Character(len=10), dimension(3)  :: nsymb
4134
       Character(len=3),  dimension(3)  :: ssymb
4135
       real(kind=cp),     parameter     :: epss=0.001
4136

    
4137
       if(present(ord) .and. present(ss) .and. present(att)) then
4138
         order=ord
4139
         ss_ptr(1:order) = ss(1:ord)
4140
         atr(:,1:order)  = att(:,1:ord)
4141
       else
4142
         call get_stabilizer(x,Spgr,order,ss_ptr,atr)
4143
       end if
4144

    
4145
       !If codes were not assigned with explicit number
4146
       !attribute numbers bigger than initial Codini
4147

    
4148
       do j=1,3
4149
          if(Icodes(j) < 1  .and. abs(multip(j)) > epss)  then
4150
               codini = codini+1
4151
               Icodes(j) = codini
4152
          end if
4153
       end do
4154

    
4155
       ssymb=(/"  x","  y","  z"/)
4156

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

    
4159
       if(order > 1 ) then  !A constraint in position must exist
4160

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

    
4163
          do k=2,order
4164
           symbol=" "
4165
           Rsym=Spgr%SymOp(ss_ptr(k))%Rot
4166
           tr=Spgr%SymOp(ss_ptr(k))%tr + atr(:,k)
4167
           call Get_SymSymb(Rsym,tr,Sym_Symb)
4168
           call symmetry_symbol(Sym_Symb,tsymbol)
4169
             i=index(tsymbol,";")
4170
             if(i /= 0) then
4171
               symbol=tsymbol(1:i-1)
4172
               call Read_Xsym(tsymbol(i+1:),1,Rsym,Tr,.false.)
4173
               if(sum(abs(x-tr)) < epss) then
4174
                  ssymb=(/"  0","  0","  0"/)
4175
                  if(present(Ipr)) then
4176
                    write(unit=Ipr,fmt="(a,i2,a,t20,a,t55,a,t90,4a)") "     Operator ",k,": ", &
4177
                    trim(Sym_Symb),trim(tsymbol),"  ssymb:" ,(ssymb(j)//"  ",j=1,3)
4178
                  end if
4179
                  cycle
4180
               end if
4181
             else
4182
               symbol=tsymbol
4183
             end if
4184
             ipar=index(symbol,")")              !Translation element appears before position
4185
             L =index(symbol(ipar+1:)," ")+ipar  !Position of the first blank after translation
4186
             L1=index(symbol(ipar+1:),",")+ipar  !Position of the first comma after translation
4187
             L2=index(symbol(L1+1:),",")+L1      !Position of the second comma
4188
             if(L1 == 0) L1=1
4189
             if(L2 == 0) L2=1
4190
             if(L  == 0) L=1
4191

    
4192
             !Construct a new symbol that estabish automatically the constraints
4193
             nsymb = (/symbol(L+1:L1-1),symbol(L1+1:L2-1),symbol(L2+1:)/)
4194

    
4195
             do i=1,3
4196
                 do j=1,10  !Delete unwanted symbols (keep only x,y,z,2 and -
4197
                    if(nsymb(i)(j:j) == " ") cycle
4198
                    if(nsymb(i)(j:j) /= "x" .and. nsymb(i)(j:j) /= "y" .and. &
4199
                       nsymb(i)(j:j) /= "z" .and. nsymb(i)(j:j) /= "-" .and. &
4200
                       nsymb(i)(j:j) /= "2" ) nsymb(i)(j:j)=" "
4201
                 end do
4202
                 if(len_trim(nsymb(i))  == 0 .or. (index(nsymb(i),"x") == 0 .and. &
4203
                    index(nsymb(i),"y") == 0 .and. index(nsymb(i),"z") == 0  ) ) then
4204
                    ssymb(i)="  0"
4205
                    cycle
4206
                 end if
4207
                 !Now remove 2s on the right of x,y, or z
4208
                 j1=index(nsymb(i),"2")
4209
                 if( j1 /= 0) then
4210
                    if(len_trim(nsymb(i)) == j1) nsymb(i)=nsymb(i)(1:j1-1)
4211
                 end if
4212
                 !Now remove -s on the right of x,y, or z
4213
                 j1=index(nsymb(i),"-")
4214
                 if( j1 /= 0) then
4215
                    if(len_trim(nsymb(i)) == j1) nsymb(i)=nsymb(i)(1:j1-1)
4216
                 end if
4217
                 nsymb(i)= adjustl(nsymb(i))
4218
             end do
4219

    
4220
             if(ssymb(1) /= "  0" .and. ssymb(1) /= "  a") then
4221
                ssymb(1)= nsymb(1)
4222
                ssymb(1)= adjustr(ssymb(1))
4223
             end if
4224

    
4225
             if(ssymb(2) /= "  0" .and. ssymb(2) /= "  a" .and. ssymb(2) /= "  b" .and. &
4226
                ssymb(2) /= " -a" .and. ssymb(2) /= " 2a"   ) then
4227
                ssymb(2) = nsymb(2)
4228
                ssymb(2) = adjustr(ssymb(2))
4229
             end if
4230

    
4231
             if(ssymb(3) /= "  0" .and. ssymb(3) /= "  a" .and. ssymb(3) /= "  b" .and. &
4232
                ssymb(3) /= "  c" .and. ssymb(3) /= " 2a" .and. ssymb(3) /= " 2b" .and. &
4233
                ssymb(3) /= " -a" .and. ssymb(3) /= " -b") then
4234
                ssymb(3) = nsymb(3)
4235
                ssymb(3) = adjustr(ssymb(3))
4236
             end if
4237

    
4238
             do i=1,3
4239
                if(ssymb(i)(3:3) == "x")  ssymb(i)(3:3) = "a"
4240
             end do
4241
             do i=1,3
4242
                if(ssymb(i)(3:3) == "y")  ssymb(i)(3:3) = "b"
4243
             end do
4244
             do i=1,3
4245
                if(ssymb(i)(3:3) == "z")  ssymb(i)(3:3) = "c"
4246
             end do
4247
             if(present(Ipr)) then
4248
                write(unit=Ipr,fmt="(a,i2,a,t20,a,t55,a,t90,4a)") "     Operator ",k,": ", &
4249
                trim(Sym_Symb),trim(tsymbol),"  Ssymb:" ,(ssymb(j)//"  ",j=1,3)
4250
             end if
4251

    
4252
          end do !do k=1,order  over operators of the stabilizer
4253

    
4254
       else
4255
         ssymb=(/"  a","  b","  c"/)
4256

    
4257
       end if  !order > 1
4258

    
4259
       do i=1,3                  !Fixing codes
4260
         if(ssymb(i)=="  0") then
4261
           Icodes(i)=0
4262
           multip(i)=0.0
4263
         end if
4264
       end do
4265

    
4266
       if(index(ssymb(1),"a") /= 0) then
4267

    
4268
         do i=2,3  !Fixing codes
4269
           if(index(ssymb(i),"-a") /= 0) then
4270
             Icodes(i)=Icodes(1)
4271
             multip(i)=-multip(1)
4272
           else if(index(ssymb(i),"a") /= 0) then
4273
             Icodes(i)=Icodes(1)
4274
             multip(i)=multip(1)
4275

    
4276
             if(index(ssymb(i),"2") /= 0) then
4277
               multip(i)=2.0* multip(1)
4278
             else if(index(ssymb(1),"2") /= 0) then
4279
               multip(i)=0.5* multip(1)
4280
             end if
4281

    
4282
           end if
4283
         end do
4284
       else  !the x-coordinate is fixed, analyse y and z
4285
         if(index(ssymb(2),"b") /= 0 .and. index(ssymb(3),"b") /= 0) then
4286
           Icodes(3)=Icodes(2)
4287
           if(ssymb(2) == ssymb(3)) then
4288
             multip(3)= multip(2)
4289
           else if(ssymb(3) == " -b" .and. ssymb(2) == "  b") then
4290
             multip(3)= -multip(2)
4291
           else if(ssymb(3) == "  b" .and. ssymb(2) == " -b") then
4292
             multip(3)= -multip(2)
4293
           end if
4294
         end if
4295

    
4296
       end if !if(index(ssymb(1),"a") /= 0)
4297

    
4298
       do j=1,3
4299
         if(abs(multip(j)) < epss) then
4300
           Icodes(j) = 0
4301
         end if
4302
       end do
4303
       if(present(Ipr)) then
4304
         write(unit=Ipr,fmt="(a,3i5)")    "     Codes positions: ",Icodes
4305
         write(unit=Ipr,fmt="(a,3f5.1)")  "     Multipliers    : ",multip
4306
         write(unit=Ipr,fmt="(5a)")       "     Codes   string : ( ",(ssymb(j),j=1,3) ," )"
4307
       end if
4308
       return
4309
    End Subroutine Get_Atompos_Ctr
4310

    
4311
    !!--++
4312
    !!--++ Subroutine Get_ConCodes_Line(Line,FAtom/MolCrys/Molec)
4313
    !!--++    character(len=*),             intent(in)     :: Line
4314
    !!--++    integer,                      intent(in)     :: Nat
4315
    !!--++    type(Atom_List_Type),         intent(in out) :: FAtom
4316
    !!--++    or
4317
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
4318
    !!--++    or
4319
    !!--++    type(molecule_type),          intent(in out) :: Molec
4320
    !!--++
4321
    !!--++    (Private)
4322
    !!--++    Get the Constraints relations
4323
    !!--++
4324
    !!--++ Update: March - 2005
4325
    !!
4326

    
4327
    !!--++
4328
    !!--++ Subroutine Get_ConCodes_Line_FAtom(Line,FAtom)
4329
    !!--++    character(len=*),         intent(in)     :: Line
4330
    !!--++    integer,                  intent(in)     :: Nat
4331
    !!--++    type(Atom_List_Type),     intent(in out) :: FAtom
4332
    !!--++
4333
    !!--++ Overloaded
4334
    !!--++ Get the Constraints relations
4335
    !!--++
4336
    !!--++ Update: March - 2005
4337
    !!
4338
    Subroutine Get_ConCodes_Line_FAtom(Line,FAtom)
4339
       !---- Arguments ----!
4340
       character(len=*),     intent(in)     :: Line
4341
       type(Atom_List_Type), intent(in out) :: FAtom
4342

    
4343
       !---- Local variables ----!
4344
       character(len=20), dimension(30) :: label
4345
       integer                          :: j,ic,n,na,nb,nc,nd,npos
4346
       integer                          :: nl,nl2,iv
4347
       integer, dimension(1)            :: ivet
4348
       real(kind=cp)                    :: fac_0,fac_1
4349
       real(kind=cp),dimension(1)       :: vet
4350

    
4351
       call init_err_refcodes()
4352

    
4353
       nl=0
4354
       call getword(line,label,ic)
4355
       if (ic < 2) then
4356
          err_refcodes=.true.
4357
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
4358
          return
4359
       end if
4360

    
4361
       !---- Set Futher Information----!
4362
       !---- Na is the number of atom on List
4363
       !---- Nb is the key (X,Y,Z,Occ,...)
4364
       !---- Fac0 is the multiplier
4365
       !---- Nl is the number of refinement parameter
4366

    
4367
       npos=index(label(1),"_")
4368
       if (npos ==0) then
4369
          err_refcodes=.true.
4370
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4371
          return
4372
       end if
4373

    
4374
       na=0
4375
       do j=1,FAtom%Natoms
4376
          if (u_case(FAtom%atom(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4377
             na=j
4378
             exit
4379
          end if
4380
       end do
4381
       if (na == 0) then
4382
          err_refcodes=.true.
4383
          ERR_RefCodes_Mess=" Atom label not found for "//trim(line)
4384
          return
4385
       end if
4386

    
4387
       nb=0
4388
       do j=1,ncode
4389
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
4390
             nb=j
4391
             exit
4392
          end if
4393
       end do
4394
       if (nb == 0) then
4395
          err_refcodes=.true.
4396
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
4397
          return
4398
       end if
4399

    
4400
       select case (nb)
4401
          case ( 1:3)
4402
             fac_0=FAtom%atom(na)%mx(nb)
4403
                nl=FAtom%atom(na)%lx(nb)
4404
          case ( 4)
4405
             fac_0=FAtom%atom(na)%mbiso
4406
                nl=FAtom%atom(na)%lbiso
4407
          case ( 5)
4408
             fac_0=FAtom%atom(na)%mocc
4409
                nl=FAtom%atom(na)%locc
4410
          case ( 6:11)
4411
             fac_0=FAtom%atom(na)%mu(nb-5)
4412
                nl=FAtom%atom(na)%lu(nb-5)
4413
          case (12)
4414
             fac_0=FAtom%atom(na)%mu(1)
4415
          case (13:)
4416
             err_refcodes=.true.
4417
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4418
             return
4419
       end select ! nb
4420

    
4421
       if (nb < ncode) then
4422
          if (nl == 0) then
4423
             err_refcodes=.true.
4424
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
4425
             return
4426
          end if
4427
       end if
4428

    
4429
       !---- Set the rest elements in Contsraints ----!
4430
       n=1
4431
       do
4432
          n=n+1
4433
          if (n > ic) exit
4434

    
4435
          npos=index(label(n),"_")
4436
          if (npos ==0) then
4437
             err_refcodes=.true.
4438
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
4439
             return
4440
          end if
4441

    
4442
          nc=0
4443
          do j=1,FAtom%Natoms
4444
             if (u_case(FAtom%atom(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
4445
                nc=j
4446
                exit
4447
             end if
4448
          end do
4449
          if (nc == 0) then
4450
             err_refcodes=.true.
4451
             ERR_RefCodes_Mess=" Atom label not found for "//trim(label(n))
4452
             return
4453
          end if
4454

    
4455
          nd=0
4456
          do j=1,ncode
4457
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
4458
                nd=j
4459
                exit
4460
             end if
4461
          end do
4462
          if (nd == 0) then
4463
             err_refcodes=.true.
4464
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
4465
             return
4466
          end if
4467

    
4468
          !---- Is there a new multiplier?
4469
          n=n+1
4470
          call getnum(label(n),vet,ivet,iv)
4471
          if (iv == 1) then
4472
             fac_1=vet(1)
4473
          else
4474
             fac_1=fac_0
4475
             n=n-1
4476
          end if
4477

    
4478
          select case (nd)
4479
             case ( 1:3)
4480
                nl2=FAtom%atom(nc)%lx(nd)
4481
                call Delete_refCodes(nl2,FAtom)
4482
                FAtom%atom(nc)%mx(nd)=fac_1
4483
                FAtom%atom(nc)%lx(nd)=nl
4484
             case ( 4)
4485
                nl2=FAtom%atom(nc)%lbiso
4486
                call Delete_refCodes(nl2,FAtom)
4487
                FAtom%atom(nc)%mbiso=fac_1
4488
                FAtom%atom(nc)%lbiso=nl
4489
             case ( 5)
4490
                nl2=FAtom%atom(nc)%locc
4491
                call Delete_refCodes(nl2,FAtom)
4492
                FAtom%atom(nc)%mocc=fac_1
4493
                FAtom%atom(nc)%locc=nl
4494
             case ( 6:11)
4495
                nl2=FAtom%atom(nc)%lu(nd-5)
4496
                call Delete_refCodes(nl2,FAtom)
4497
                FAtom%atom(nc)%mu(nd-5)=fac_1
4498
                FAtom%atom(nc)%lu(nd-5)=nl
4499
             case (12)
4500
                do j=1,6
4501
                   nl2=FAtom%atom(nc)%lu(j)
4502
                   call Delete_refCodes(nl2,FAtom)
4503
                   FAtom%atom(nc)%mu(j)=fac_1
4504
                   FAtom%atom(nc)%lu(j)=FAtom%atom(na)%lu(j)
4505
                end do
4506
                np_cons=np_cons+5
4507
             case (13:)
4508
                err_refcodes=.true.
4509
                ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4510
                return
4511
          end select ! nb
4512

    
4513
          np_cons=np_cons+1
4514

    
4515
       end do
4516

    
4517
       return
4518
    End Subroutine Get_ConCodes_Line_FAtom
4519

    
4520
    !!--++
4521
    !!--++ Subroutine Get_ConCodes_Line_FmAtom(Line,FmAtom)
4522
    !!--++    character(len=*),         intent(in)     :: Line
4523
    !!--++    integer,                  intent(in)     :: Nat
4524
    !!--++    type(mAtom_List_Type),    intent(in out) :: FmAtom
4525
    !!--++
4526
    !!--++ Get the Magnetic Constraints Relations: Presently only 'equal'
4527
    !!--++ mag clone of Get_ConCodes_Line_FAtom
4528
    !!--++ Created: December - 2011
4529
    !!
4530
    Subroutine Get_ConCodes_Line_FmAtom(Line,FmAtom)
4531
       !---- Arguments ----!
4532
       character(len=*),     intent(in)     :: Line
4533
       type(mAtom_List_Type),intent(in out) :: FmAtom
4534

    
4535
       !---- Local variables ----!
4536
       character(len=20), dimension(30) :: label
4537
       integer                          :: j,ic,n,na,nb,nc,nd,npos
4538
       integer                          :: nl,nl2,iv,ik
4539
       integer, dimension(1)            :: ivet
4540
       real(kind=cp)                    :: fac_0,fac_1
4541
       real(kind=cp),dimension(1)       :: vet
4542

    
4543
       call init_err_refcodes()
4544

    
4545
       nl=0
4546
       call getword(line,label,ic)
4547
       if (ic < 2) then
4548
          err_refcodes=.true.
4549
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
4550
          return
4551
       end if
4552

    
4553
       !---- Set Futher Information----!
4554
       !---- Na is the number of atom on List
4555
       !---- Nb is the number of keys (Rx,Ry,Rz,Ix,Iy,Iz,MagPh...)
4556
       !---- Fac0 is the multiplier
4557
       !---- Nl is the number of refinement parameter
4558

    
4559
       npos=index(label(1),"_")
4560
       if (npos ==0) then
4561
          err_refcodes=.true.
4562
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4563
          return
4564
       end if
4565

    
4566
       na=0
4567
       do j=1,FmAtom%Natoms
4568
          if (u_case(FmAtom%atom(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4569
             na=j
4570
             exit
4571
          end if
4572
       end do
4573
       if (na == 0) then
4574
          err_refcodes=.true.
4575
          ERR_RefCodes_Mess=" Atom label not found for "//trim(line)
4576
          return
4577
       end if
4578

    
4579
       nb=0
4580
       do j=1,ncode
4581
          if (u_case(label(1)(1:npos))==u_case(trim(mcode_nam(j)))) then
4582
             nb=j
4583
             exit
4584
          end if
4585
       end do
4586
       if (nb == 0) then
4587
          err_refcodes=.true.
4588
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
4589
          return
4590
       end if
4591

    
4592
       !---- Get im, number of the magnetic matrices/irrep set
4593
       ik=FmAtom%atom(na)%nvk
4594
       
4595
       select case (nb)
4596
          case ( 1:3)
4597
             fac_0=FmAtom%atom(na)%mSkR(nb,ik)
4598
                nl=FmAtom%atom(na)%lSkR(nb,ik)
4599
          case ( 4:6)
4600
             fac_0=FmAtom%atom(na)%mSkI(nb-3,ik)
4601
                nl=FmAtom%atom(na)%lSkI(nb-3,ik)
4602
          case ( 7:9)
4603
             fac_0=FmAtom%atom(na)%mSkR(nb-6,ik)
4604
                nl=FmAtom%atom(na)%lSkR(nb-6,ik)
4605
          case (10:12)
4606
             fac_0=FmAtom%atom(na)%mSkI(nb-9,ik)
4607
                nl=FmAtom%atom(na)%lSkI(nb-9,ik)
4608
          case (13)
4609
             fac_0=FmAtom%atom(na)%mmphas(ik)
4610
                nl=FmAtom%atom(na)%lmphas(ik)
4611
           case (14:22)
4612
             fac_0=FmAtom%atom(na)%mbas(nb-13,ik)
4613
                nl=FmAtom%atom(na)%lbas(nb-13,ik)
4614
         case (23:)
4615
             err_refcodes=.true.
4616
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4617
             return
4618
       end select ! nb
4619

    
4620
       if (nb < ncode) then
4621
          if (nl == 0) then
4622
             err_refcodes=.true.
4623
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
4624
             return
4625
          end if
4626
       end if
4627

    
4628
       !---- Set the rest elements in Contsraints ----!
4629
       n=1
4630
       do
4631
          n=n+1
4632
          if (n > ic) exit
4633

    
4634
          npos=index(label(n),"_")
4635
          if (npos ==0) then
4636
             err_refcodes=.true.
4637
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
4638
             return
4639
          end if
4640

    
4641
          nc=0
4642
          do j=1,FmAtom%Natoms
4643
             if (u_case(FmAtom%atom(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
4644
                nc=j
4645
                exit
4646
             end if
4647
          end do
4648
          if (nc == 0) then
4649
             err_refcodes=.true.
4650
             ERR_RefCodes_Mess=" Atom label not found for "//trim(label(n))
4651
             return
4652
          end if
4653

    
4654
          nd=0
4655
          do j=1,ncode
4656
             if (u_case(label(n)(1:npos))==u_case(trim(mcode_nam(j)))) then
4657
                nd=j
4658
                exit
4659
             end if
4660
          end do
4661
          if (nd == 0) then
4662
             err_refcodes=.true.
4663
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
4664
             return
4665
          end if
4666

    
4667
          !---- Is there a new multiplier?
4668
          n=n+1
4669
          call getnum(label(n),vet,ivet,iv)
4670
          if (iv == 1) then
4671
             fac_1=vet(1)
4672
          else
4673
             fac_1=fac_0
4674
             n=n-1
4675
          end if
4676

    
4677
          select case (nd)
4678
             case ( 1:3)
4679
                nl2=FmAtom%atom(nc)%lSkR(nd,ik)
4680
                call Delete_refCodes(nl2,FmAtom)
4681
                FmAtom%atom(nc)%mSkR(nd,ik)=fac_1
4682
                FmAtom%atom(nc)%lSkR(nd,ik)=nl
4683
             case ( 4:6)
4684
                nl2=FmAtom%atom(nc)%lSkI(nd-3,ik)
4685
                call Delete_refCodes(nl2,FmAtom)
4686
                FmAtom%atom(nc)%mSkI(nd-3,ik)=fac_1
4687
                FmAtom%atom(nc)%lSkI(nd-3,ik)=nl
4688
             case ( 7:9)
4689
                nl2=FmAtom%atom(nc)%lSkR(nd-6,ik)
4690
                call Delete_refCodes(nl2,FmAtom)
4691
                FmAtom%atom(nc)%mSkR(nd-6,ik)=fac_1
4692
                FmAtom%atom(nc)%lSkR(nd-6,ik)=nl
4693
             case (10:12)
4694
                nl2=FmAtom%atom(nc)%lSkI(nd-9,ik)
4695
                call Delete_refCodes(nl2,FmAtom)
4696
                FmAtom%atom(nc)%mSkI(nd-9,ik)=fac_1
4697
                FmAtom%atom(nc)%lSkI(nd-9,ik)=nl
4698
             case (13)
4699
                nl2=FmAtom%atom(nc)%lmphas(ik)
4700
                call Delete_refCodes(nl2,FmAtom)
4701
                FmAtom%atom(nc)%mmphas(ik)=fac_1
4702
                FmAtom%atom(nc)%lmphas(ik)=nl
4703
             case (14:22)
4704
                nl2=FmAtom%atom(nc)%lbas(nd-13,ik)
4705
                call Delete_refCodes(nl2,FmAtom)
4706
                FmAtom%atom(nc)%mbas(nd-13,ik)=fac_1
4707
                FmAtom%atom(nc)%lbas(nd-13,ik)=nl
4708
             case (23:)
4709
                err_refcodes=.true.
4710
                ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
4711
                return
4712
          end select ! nb
4713

    
4714
          np_cons=np_cons+1
4715

    
4716
       end do
4717

    
4718
       return
4719
    End Subroutine Get_ConCodes_Line_FmAtom
4720

    
4721
    !!--++
4722
    !!--++ Subroutine Get_ConCodes_Line_Molcrys(Line,Molcrys)
4723
    !!--++    character(len=*),             intent(in)     :: Line
4724
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
4725
    !!--++
4726
    !!--++ Overloaded
4727
    !!--++ Get the Constraints relations
4728
    !!--++
4729
    !!--++ Update: March - 2005
4730
    !!
4731
    Subroutine Get_ConCodes_Line_Molcrys(Line,Molcrys)
4732
       !---- Arguments ----!
4733
       character(len=*),             intent(in)     :: Line
4734
       type(molecular_Crystal_type), intent(in out) :: MolCrys
4735

    
4736
       !---- Loval variables ----!
4737
       character(len=5)                 :: car
4738
       character(len=20), dimension(30) :: label
4739
       integer                          :: i,j,ic,n,na,naa,nb,nc,ncc,nd
4740
       integer                          :: npos, nposm, nmol1,nmol2
4741
       integer                          :: nl,nl2,iv
4742
       integer, dimension(1)            :: ivet
4743
       real(kind=cp)                    :: fac_0,fac_1
4744
       real(kind=cp),dimension(1)       :: vet
4745

    
4746
       call init_err_refcodes()
4747

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

    
4755
       !---- Set Father Information ----!
4756
       !---- Na is the number of atom on List
4757
       !---- Nb is the key (X,Y,Z,Occ,...)
4758
       !---- Fac0 is the multiplier
4759
       !---- Nl is the number of refinement parameter
4760
       npos=index(label(1),"_")
4761
       if (npos ==0) then
4762
          err_refcodes=.true.
4763
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
4764
          return
4765
       end if
4766
       nposm=index(u_case(label(1)),"MOL")
4767
       if (nposm /= 0) then
4768
          car=adjustl(label(1)(nposm+3:))
4769
          if (car(1:2) == "  ") then
4770
             nmol1=0
4771
          else
4772
             read(unit=car,fmt="(i2)") nmol1
4773
          end if
4774
       else
4775
          nmol1=-1
4776
       end if
4777

    
4778
       na=0
4779
       do j=1,molcrys%n_free
4780
          if (u_case(molcrys%atm(j)%lab) == u_case(label(1)(npos+1:npos+6))) then
4781
             na=j
4782
             exit
4783
          end if
4784
       end do
4785
       if (na == 0) then
4786
          do i=1,molcrys%n_mol
4787
             do j=1,molcrys%mol(i)%natoms
4788
                if (u_case(molcrys%mol(i)%Atname(j)) == u_case(label(1)(npos+1:npos+6))) then
4789
                   if (j > 1) then
4790
                      na=molcrys%n_free+sum(molcrys%mol(1:j-1)%natoms)+j
4791
                   else
4792
                      na=molcrys%n_free+j
4793
                   end if
4794
                   exit
4795
                end if
4796
             end do
4797
          end do
4798
       end if
4799

    
4800
       nb=0
4801
       do j=1,ncode
4802
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
4803
             nb=j
4804
             exit
4805
          end if
4806
       end do
4807
       if (nb == 0) then
4808
          err_refcodes=.true.
4809
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
4810
          return
4811
       end if
4812

    
4813
       !---- Checking ----!
4814
       if (nb <= 12 .and. na==0) then
4815
          err_refcodes=.true.
4816
          ERR_RefCodes_Mess="Incompatible option for parameter name: "//trim(label(1))
4817
          return
4818
       end if
4819

    
4820
       nl=0
4821
       select case (nb)
4822
          case ( 1:3)
4823
             !---- X_, Y_, Z_ ----!
4824
             if (na <= molcrys%n_free) then
4825
                fac_0=molcrys%atm(na)%mx(nb)
4826
                   nl=molcrys%atm(na)%lx(nb)
4827
             else
4828
                naa=na-molcrys%n_free
4829
                do i=1,molcrys%n_mol
4830
                   if (naa > molcrys%mol(i)%natoms) then
4831
                      naa=naa-molcrys%mol(i)%natoms
4832
                      cycle
4833
                   end if
4834
                   fac_0=molcrys%mol(i)%mI_coor(nb,naa)
4835
                      nl=molcrys%mol(i)%lI_coor(nb,naa)
4836
                end do
4837
             end if
4838

    
4839
          case ( 4)
4840
             !---- Biso_ ----!
4841
             if (na <= molcrys%n_free) then
4842
                fac_0=molcrys%atm(na)%mbiso
4843
                   nl=molcrys%atm(na)%lbiso
4844
             else
4845
                naa=na-molcrys%n_free
4846
                do i=1,molcrys%n_mol
4847
                   if (naa > molcrys%mol(i)%natoms) then
4848
                      naa=naa-molcrys%mol(i)%natoms
4849
                      cycle
4850
                   end if
4851
                   fac_0=molcrys%mol(i)%mbiso(naa)
4852
                      nl=molcrys%mol(i)%lbiso(naa)
4853
                end do
4854
             end if
4855

    
4856
          case ( 5)
4857
             !---- Occ_ ----!
4858
             if (na <= molcrys%n_free) then
4859
                fac_0=molcrys%atm(na)%mocc
4860
                   nl=molcrys%atm(na)%locc
4861
             else
4862
                naa=na-molcrys%n_free
4863
                do i=1,molcrys%n_mol
4864
                   if (naa > molcrys%mol(i)%natoms) then
4865
                      naa=naa-molcrys%mol(i)%natoms
4866
                      cycle
4867
                   end if
4868
                   fac_0=molcrys%mol(i)%mocc(naa)
4869
                      nl=molcrys%mol(i)%locc(naa)
4870
                end do
4871
             end if
4872

    
4873
          case ( 6:11)
4874
             !---- B11_, ..., B23_ ----!
4875
             if (na <= molcrys%n_free) then
4876
                fac_0=molcrys%atm(na)%mu(nb-5)
4877
                   nl=molcrys%atm(na)%lu(nb-5)
4878
             else
4879
                err_refcodes=.true.
4880
                ERR_RefCodes_Mess="Option no valid"
4881
                return
4882
             end if
4883

    
4884
          case (12)
4885
             !---- Banis_ ----!
4886
             if (na <= molcrys%n_free) then
4887
                fac_0=molcrys%atm(na)%mu(1)
4888
             else
4889
                err_refcodes=.true.
4890
                ERR_RefCodes_Mess="Option no valid"
4891
                return
4892
             end if
4893

    
4894
          case (13:15)
4895
             !---- Xc_, Yc_, Zc_ ----!
4896
             select case (nmol1)
4897
                case (-1)
4898
                   err_refcodes=.true.
4899
                   ERR_RefCodes_Mess="Option no valid"
4900
                   return
4901

    
4902
                case (0)
4903
                   fac_0=molcrys%mol(1)%mxcentre(nb-12)
4904
                   nl=molcrys%mol(1)%lxcentre(nb-12)
4905

    
4906
                case (1:)
4907
                   fac_0=molcrys%mol(nmol1)%mxcentre(nb-12)
4908
                   nl=molcrys%mol(nmol1)%lxcentre(nb-12)
4909
             end select
4910

    
4911
          case (16:18)
4912
             !---- Theta_ , Phi_, Chi_ ----!
4913
             select case (nmol1)
4914
                case (-1)
4915
                   err_refcodes=.true.
4916
                   ERR_RefCodes_Mess="Option no valid"
4917
                   return
4918

    
4919
                case (0)
4920
                   fac_0=molcrys%mol(1)%mOrient(nb-15)
4921
                   nl=molcrys%mol(1)%lOrient(nb-15)
4922

    
4923
                case (1)
4924
                   fac_0=molcrys%mol(nmol1)%mOrient(nb-15)
4925
                   nl=molcrys%mol(nmol1)%lOrient(nb-15)
4926
             end select
4927

    
4928
          case (19:21)
4929
             !!! Not yet Implemented !!!
4930
       end select ! nb
4931

    
4932
       if (nb < ncode) then
4933
          if (nl == 0) then
4934
             err_refcodes=.true.
4935
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
4936
             return
4937
          end if
4938
       end if
4939

    
4940
       !---- Set the rest elements in Constraints ----!
4941
       n=1
4942
       do
4943
          n=n+1
4944
          if (n > ic) exit
4945

    
4946
          npos=index(label(n),"_")
4947
          if (npos ==0) then
4948
             err_refcodes=.true.
4949
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
4950
             return
4951
          end if
4952
          nposm=index(u_case(label(n)),"MOL")
4953
          if (nposm /= 0) then
4954
             car=adjustl(label(n)(nposm+3:))
4955
             if (car(1:2) == "  ") then
4956
                nmol2=0
4957
             else
4958
                read(unit=car,fmt="(i2)") nmol2
4959
             end if
4960
          else
4961
             nmol2=-1
4962
          end if
4963

    
4964
          nc=0
4965
          do j=1,molcrys%n_free
4966
             if (u_case(molcrys%atm(j)%lab) == u_case(label(n)(npos+1:npos+6))) then
4967
                nc=j
4968
                exit
4969
             end if
4970
          end do
4971
          if (nc == 0) then
4972
             do i=1,molcrys%n_mol
4973
                do j=1,molcrys%mol(i)%natoms
4974
                   if (u_case(molcrys%mol(i)%Atname(j)) == u_case(label(n)(npos+1:npos+6))) then
4975
                      if (j > 1) then
4976
                         nc=molcrys%n_free+sum(molcrys%mol(1:j-1)%natoms)+j
4977
                      else
4978
                         nc=molcrys%n_free+j
4979
                      end if
4980
                      exit
4981
                   end if
4982
                end do
4983
             end do
4984
          end if
4985

    
4986
          nd=0
4987
          do j=1,ncode
4988
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
4989
                nd=j
4990
                exit
4991
             end if
4992
          end do
4993
          if (nd == 0) then
4994
             err_refcodes=.true.
4995
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
4996
             return
4997
          end if
4998

    
4999
          !---- Checking ----!
5000
          if (nd <= 12 .and. nc==0) then
5001
             err_refcodes=.true.
5002
             ERR_RefCodes_Mess="Incompatible option for parameter name: "//trim(label(n))
5003
             return
5004
          end if
5005

    
5006
          !---- Is there a new multiplier? ----!
5007
          n=n+1
5008
          call getnum(label(n),vet,ivet,iv)
5009
          if (iv == 1) then
5010
             fac_1=vet(1)
5011
          else
5012
             fac_1=fac_0
5013
             n=n-1
5014
          end if
5015

    
5016
          select case (nd)
5017
             case ( 1:3)
5018
                !---- X_, Y_, Z_ ----!
5019
                if (nc <= molcrys%n_free) then
5020
                   nl2=molcrys%atm(nc)%lx(nd)
5021
                   call Delete_RefCodes(nl2,molcrys)
5022
                   molcrys%atm(nc)%mx(nd)=fac_1
5023
                   molcrys%atm(nc)%lx(nd)=nl
5024
                else
5025
                   ncc=nc-molcrys%n_free
5026
                   do i=1,molcrys%n_mol
5027
                      if (ncc > molcrys%mol(i)%natoms) then
5028
                         ncc=ncc-molcrys%mol(i)%natoms
5029
                         cycle
5030
                      end if
5031
                      nl2=molcrys%mol(i)%lI_coor(nd,ncc)
5032
                      call Delete_RefCodes(nl2,molcrys)
5033
                      molcrys%mol(i)%mI_coor(nd,ncc)=fac_1
5034
                      molcrys%mol(i)%lI_coor(nd,ncc)=nl
5035
                   end do
5036
                end if
5037

    
5038
             case ( 4)
5039
                !---- Biso_ ----!
5040
                if (nc <= molcrys%n_free) then
5041
                   nl2=molcrys%atm(nc)%lbiso
5042
                   call Delete_RefCodes(nl2,molcrys)
5043
                   molcrys%atm(nc)%mbiso=fac_1
5044
                   molcrys%atm(nc)%lbiso=nl
5045
                else
5046
                   ncc=nc-molcrys%n_free
5047
                   do i=1,molcrys%n_mol
5048
                      if (ncc > molcrys%mol(i)%natoms) then
5049
                         ncc=ncc-molcrys%mol(i)%natoms
5050
                         cycle
5051
                      end if
5052
                      nl2=molcrys%mol(i)%lbiso(ncc)
5053
                      call Delete_RefCodes(nl2,molcrys)
5054
                      molcrys%mol(i)%mbiso(ncc)=fac_1
5055
                      molcrys%mol(i)%lbiso(ncc)=nl
5056
                   end do
5057
                end if
5058

    
5059
             case ( 5)
5060
                !---- Occ_ ----!
5061
                if (nc <= molcrys%n_free) then
5062
                   nl2=molcrys%atm(nc)%locc
5063
                   call Delete_RefCodes(nl2,molcrys)
5064
                   molcrys%atm(nc)%mocc=fac_1
5065
                   molcrys%atm(nc)%locc=nl
5066
                else
5067
                   ncc=nc-molcrys%n_free
5068
                   do i=1,molcrys%n_mol
5069
                      if (ncc > molcrys%mol(i)%natoms) then
5070
                         ncc=ncc-molcrys%mol(i)%natoms
5071
                         cycle
5072
                      end if
5073
                      nl2=molcrys%mol(i)%locc(ncc)
5074
                      call Delete_RefCodes(nl2,molcrys)
5075
                      molcrys%mol(i)%mocc(ncc)=fac_1
5076
                      molcrys%mol(i)%locc(ncc)=nl
5077
                   end do
5078
                end if
5079

    
5080
             case ( 6:11)
5081
                !---- B11_, ...., B23_ ----!
5082
                if (nc <= molcrys%n_free) then
5083
                   nl2=molcrys%atm(nc)%lu(nd-5)
5084
                   call Delete_RefCodes(nl2,molcrys)
5085
                   molcrys%atm(nc)%mu(nd-5)=fac_1
5086
                   molcrys%atm(nc)%lu(nd-5)=nl
5087
                else
5088
                   err_refcodes=.true.
5089
                   ERR_RefCodes_Mess="Option no valid"
5090
                   return
5091
                end if
5092

    
5093
             case (12)
5094
                !---- Banis_ ----!
5095
                if (nc <= molcrys%n_free) then
5096
                   do j=1,6
5097
                      nl2=molcrys%atm(nc)%lu(j)
5098
                      call Delete_RefCodes(nl2,molcrys)
5099
                      molcrys%atm(nc)%mu(j)=fac_1
5100
                      molcrys%atm(nc)%lu(j)=molcrys%atm(na)%lu(j)
5101
                   end do
5102
                   np_cons=np_cons+5
5103
                else
5104
                   err_refcodes=.true.
5105
                   ERR_RefCodes_Mess="Option no valid"
5106
                   return
5107
                end if
5108

    
5109
             case (13:15)
5110
                select case (nmol2)
5111
                   case (-1)
5112
                      err_refcodes=.true.
5113
                      ERR_RefCodes_Mess="Option no valid"
5114
                      return
5115

    
5116
                   case (0)
5117
                      do i=1,molcrys%n_mol
5118
                         nl2=molcrys%mol(i)%lxcentre(nc-12)
5119
                         call Delete_RefCodes(nl2,molcrys)
5120
                         molcrys%mol(i)%mxcentre(nc-12)=fac_1
5121
                         molcrys%mol(i)%lxcentre(nc-12)=nl
5122
                      end do
5123
                      np_cons=np_cons+(i-1)
5124

    
5125
                   case (1:)
5126
                      nl2=molcrys%mol(nmol2)%lxcentre(nc-12)
5127
                      call Delete_RefCodes(nl2,molcrys)
5128
                      molcrys%mol(nmol2)%mxcentre(nc-12)=fac_1
5129
                      molcrys%mol(nmol2)%lxcentre(nc-12)=nl
5130
                end select
5131

    
5132
             case (16:18)
5133
                select case (nmol2)
5134
                   case (-1)
5135
                      err_refcodes=.true.
5136
                      ERR_RefCodes_Mess="Option no valid"
5137
                      return
5138

    
5139
                   case (0)
5140
                      do i=1,molcrys%n_mol
5141
                         nl2=molcrys%mol(i)%lorient(nc-15)
5142
                         call Delete_RefCodes(nl2,molcrys)
5143
                         molcrys%mol(i)%morient(nc-15)=fac_1
5144
                         molcrys%mol(i)%lorient(nc-15)=nl
5145
                      end do
5146
                      np_cons=np_cons+(i-1)
5147

    
5148
                   case (1:)
5149
                      nl2=molcrys%mol(nmol2)%lorient(nc-15)
5150
                      call Delete_RefCodes(nl2,molcrys)
5151
                      molcrys%mol(nmol2)%morient(nc-15)=fac_1
5152
                      molcrys%mol(nmol2)%lorient(nc-15)=nl
5153
                end select
5154

    
5155
             case (19:21)
5156
                !!! Not yet implemented !!!
5157

    
5158
          end select ! nb
5159
          np_cons=np_cons+1
5160

    
5161
       end do
5162

    
5163
       return
5164
    End Subroutine Get_ConCodes_Line_Molcrys
5165

    
5166
    !!--++
5167
    !!--++ Subroutine Get_ConCodes_Line_Molec(Line,Molec)
5168
    !!--++    character(len=*),    intent(in)     :: Line
5169
    !!--++    type(molecule_type), intent(in out) :: Molec
5170
    !!--++
5171
    !!--++ Overloaded
5172
    !!--++ Get the Constraints relations
5173
    !!--++
5174
    !!--++ Update: March - 2005
5175
    !!
5176
    Subroutine Get_ConCodes_Line_Molec(Line,Molec)
5177
       !---- Arguments ----!
5178
       character(len=*),    intent(in)     :: Line
5179
       type(molecule_type), intent(in out) :: Molec
5180

    
5181
       !---- Loval variables ----!
5182
       character(len=20), dimension(30) :: label
5183
       integer                          :: j,ic,na,nb,nc,nd,npos!,i,naa,ncc
5184
       integer                          :: n,nl,nl2,iv
5185
       integer, dimension(1)            :: ivet
5186
       real(kind=cp)                    :: fac_0,fac_1
5187
       real(kind=cp),dimension(1)       :: vet
5188

    
5189
       call init_err_refcodes()
5190

    
5191
       call getword(line,label,ic)
5192
       if (ic < 2) then
5193
          err_refcodes=.true.
5194
          ERR_RefCodes_Mess="EQUAL keyword needs two labels: "//trim(line)
5195
          return
5196
       end if
5197

    
5198
       !---- Set Father Information ----!
5199
       !---- Na is the number of atom on List
5200
       !---- Nb is the key (X,Y,Z,Occ,...)
5201
       !---- Fac0 is the multiplier
5202
       !---- Nl is the number of refinement parameter
5203
       npos=index(label(1),"_")
5204
       if (npos ==0) then
5205
          err_refcodes=.true.
5206
          ERR_RefCodes_Mess="The name "//trim(label(1))//" does not fit any known code-name of CrysFML "
5207
          return
5208
       end if
5209

    
5210
       na=0
5211
       do j=1,molec%natoms
5212
          if (u_case(molec%Atname(j)) == u_case(label(1)(npos+1:npos+6))) then
5213
             na=j
5214
             exit
5215
          end if
5216
       end do
5217

    
5218
       nb=0
5219
       do j=1,ncode
5220
          if (u_case(label(1)(1:npos))==u_case(trim(code_nam(j)))) then
5221
             nb=j
5222
             exit
5223
          end if
5224
       end do
5225
       if (nb == 0) then
5226
          err_refcodes=.true.
5227
          ERR_RefCodes_Mess="Code-name not found for parameter name: "//trim(label(1))
5228
          return
5229
       end if
5230

    
5231
       !---- Checking ----!
5232
       if (nb < 6 .and. na == 0) then
5233
          err_refcodes=.true.
5234
          ERR_RefCodes_Mess="Incompatible relation: "//trim(label(1))
5235
          return
5236
       end if
5237

    
5238
       nl=0
5239
       select case (nb)
5240
          case ( 1:3)
5241
             !---- X_, Y_, Z_ ----!
5242
             fac_0=molec%mI_coor(nb,na)
5243
                nl=molec%lI_coor(nb,na)
5244
          case ( 4)
5245
             !---- Biso_ ----!
5246
             fac_0=molec%mbiso(na)
5247
                nl=molec%lbiso(na)
5248
          case ( 5)
5249
             !---- Occ_ ----!
5250
             fac_0=molec%mocc(na)
5251
                nl=molec%locc(na)
5252
          case ( 6:12)
5253
             !---- Anisotropic Parameters ----!
5254
             err_refcodes=.true.
5255
             ERR_RefCodes_Mess="Incompatible Code-name for parameter name: "//trim(label(1))
5256
             return
5257
          case (13:15)
5258
             !---- Xc_, Yc_, Zc_ ----!
5259
             fac_0=molec%mxcentre(nb-12)
5260
                nl=molec%lxcentre(nb-12)
5261
          case (16:18)
5262
             !---- Theta_, Phi_, Chi_ ----!
5263
             fac_0=molec%mxcentre(nb-15)
5264
                nl=molec%lxcentre(nb-15)
5265
          case (19:21)
5266
             !!! Not Yet Implemented !!!
5267
       end select ! nb
5268

    
5269
       if (nb < ncode) then
5270
          if (nl == 0) then
5271
             err_refcodes=.true.
5272
             ERR_RefCodes_Mess="No refinable parameter was selected for "//trim(label(1))
5273
             return
5274
          end if
5275
       end if
5276

    
5277
       !---- Set Others ----!
5278
       n=1
5279
       do
5280
          n=n+1
5281
          if (n > ic) exit
5282

    
5283
          npos=index(label(n),"_")
5284
          if (npos ==0) then
5285
             err_refcodes=.true.
5286
             ERR_RefCodes_Mess="No CrysFML code-name was found for "//trim(label(n))
5287
             return
5288
          end if
5289

    
5290
          nc=0
5291
          do j=1,molec%natoms
5292
             if (u_case(molec%Atname(j)) == u_case(label(n)(npos+1:npos+6))) then
5293
                nc=j
5294
                exit
5295
             end if
5296
          end do
5297

    
5298
          nd=0
5299
          do j=1,ncode
5300
             if (u_case(label(n)(1:npos))==u_case(trim(code_nam(j)))) then
5301
                nd=j
5302
                exit
5303
             end if
5304
          end do
5305
          if (nd == 0) then
5306
             err_refcodes=.true.
5307
             ERR_RefCodes_Mess="Code-name not found for "//trim(label(n))
5308
             return
5309
          end if
5310

    
5311
          n=n+1
5312
          call getnum(label(n),vet,ivet,iv)
5313
          if (iv == 1) then
5314
             fac_1=vet(1)
5315
          else
5316
             fac_1=fac_0
5317
             n=n-1
5318
          end if
5319

    
5320
          !---- Checking ----!
5321
          if (nd < 6 .and. nc == 0) then
5322
             err_refcodes=.true.
5323
             ERR_RefCodes_Mess="Incompatible relation: "//trim(label(n))
5324
             return
5325
          end if
5326

    
5327
          select case (nd)
5328
             case ( 1:3)
5329
                !---- X_, Y_, Z_ ----!
5330
                nl2=molec%lI_coor(nd,nc)
5331
                call Delete_RefCodes(nl2,molec)
5332
                molec%mI_coor(nd,nc)=fac_1
5333
                molec%lI_coor(nd,nc)=nl
5334
             case ( 4)
5335
                !---- Biso_ ----!
5336
                nl2=molec%lbiso(nc)
5337
                call Delete_RefCodes(nl2,molec)
5338
                molec%mbiso(nc)=fac_1
5339
                molec%lbiso(nc)=nl
5340
             case ( 5)
5341
                !---- Occ_ ----!
5342
                nl2=molec%locc(nc)
5343
                call Delete_RefCodes(nl2,molec)
5344
                molec%mocc(nc)=fac_1
5345
                molec%locc(nc)=nl
5346
             case ( 6:12)
5347
                err_refcodes=.true.
5348
                ERR_RefCodes_Mess="Incompatible Code-name for "//trim(label(n))
5349
                return
5350
             case (13:15)
5351
                !---- Xc_, Yc_, Zc_ ----!
5352
                nl2=molec%lxcentre(nc-12)
5353
                call Delete_RefCodes(nl2,molec)
5354
                molec%mxcentre(nc-12)=fac_1
5355
                molec%lxcentre(nc-12)=nl
5356
             case (16:18)
5357
                !---- Theta_, Phi_, Chi_ ----!
5358
                nl2=molec%lorient(nc-15)
5359
                call Delete_RefCodes(nl2,molec)
5360
                molec%morient(nc-15)=fac_1
5361
                molec%lorient(nc-15)=nl
5362
             case (19:21)
5363
                !!! not yet implemented !!!
5364
          end select ! nb
5365
          np_cons=np_cons+1
5366

    
5367
       end do
5368

    
5369
       return
5370
    End Subroutine Get_ConCodes_Line_Molec
5371

    
5372
    !!--++
5373
    !!--++ Subroutine Get_RefCodes_Line(Key,Dire,Line,FAtom/MolCrys/Molec,Spg)
5374
    !!--++    integer,                      intent(in)     :: Key
5375
    !!--++    character(len=*),             intent(in)     :: Dire
5376
    !!--++    character(len=*),             intent(in)     :: Line
5377
    !!--++    type(Atom_List_Type),         intent(in out) :: FAtom
5378
    !!--++    or
5379
    !!--++    type(molecular_Crystal_type), intent(in out) :: MolCrys
5380
    !!--++    or
5381
    !!--++    type(molecule_type),          intent(in out) :: Molec
5382
    !!--++    type(space_group_type),       intent(in)     :: Spg
5383
    !!--++
5384
    !!--++    (Private)
5385
    !!--++    Get Refinement Codes for Free atoms type
5386
    !!--++
5387
    !!--++ Update: March - 2005
5388
    !!
5389

    
5390
    !!--++
5391
    !!--++ Subroutine Get_RefCodes_Line_FAtom(Key,Dire,Line,FAtom,Spg)
5392
    !!--++    integer,                 intent(in)     :: Key
5393
    !!--++    character(len=*),        intent(in)     :: Dire
5394
    !!--++    character(len=*),        intent(in)     :: Line
5395
    !!--++    type(Atom_List_Type),    intent(in out) :: FAtom
5396
    !!--++    type(space_group_type),  intent(in)     :: Spg
5397
    !!--++
5398
    !!--++ Overloaded
5399
    !!--++ Get Refinement Codes for Free atoms type
5400
    !!--++
5401
    !!--++ Update: March - 2005
5402
    !!
5403
    Subroutine Get_RefCodes_Line_FAtom(Key,Dire,Line,FAtom,Spg)
5404
       !---- Arguments ----!
5405
       integer,                 intent(in)     :: Key
5406
       character(len=*),        intent(in)     :: Dire
5407
       character(len=*),        intent(in)     :: Line
5408
       type(Atom_List_Type),    intent(in out) :: FAtom
5409
       type(space_group_type),  intent(in)     :: Spg
5410

    
5411
       !---- Local Variables ----!
5412
       character(len=20), dimension(30) :: label
5413
       integer                          :: i,j,n,na,nb,ndir,npos,nlong,ic !,k,nc
5414
       integer                          :: icond,iv,n_ini,n_end
5415
       integer, dimension(5)            :: ivet
5416
       integer, dimension(30)           :: ilabel
5417
       real(kind=cp)                    :: x_low,x_up,x_step
5418
       real(kind=cp),dimension(5)       :: vet
5419

    
5420

    
5421
       call init_err_refcodes()
5422
       nlong=len_trim(line)
5423

    
5424
       if (nlong ==0) then
5425
          !---- Default Values ----!
5426
          do i=1,FAtom%natoms
5427
             call Fill_RefCodes(Key,Dire,i,0,0.0_cp,0.0_cp,0.0_cp,0,Fatom,Spg)
5428
          end do
5429

    
5430
       else
5431
          !---- VARY/FIX Line: [LABEL, [INF,[SUP,[STEP,[COND]]]]] ----!
5432
          ilabel=0
5433
          call getword(line,label,ic)
5434
          do i=1,ic
5435
             call getnum(label(i),vet,ivet,iv)
5436
             if (iv == 1) ilabel(i)=1
5437
          end do
5438
          ndir=count(ilabel(1:ic) < 1)
5439

    
5440
          if (ndir <=0) then
5441
             !--- [INF,[SUP,[STEP,[COND]]]] ----!
5442
             call getnum(line,vet,ivet,iv)
5443
             select case (iv)
5444
                case (1)
5445
                   x_low=vet(1)
5446
                    x_up=vet(1)
5447
                  x_step=0.0
5448
                   icond=0
5449

    
5450
                case (2)
5451
                   x_low=vet(1)
5452
                    x_up=vet(2)
5453
                  x_step=0.0
5454
                   icond=0
5455

    
5456
                case (3)
5457
                   x_low=vet(1)
5458
                    x_up=vet(2)
5459
                  x_step=vet(3)
5460
                   icond=0
5461

    
5462
                case (4)
5463
                   x_low=vet(1)
5464
                    x_up=vet(2)
5465
                  x_step=vet(3)
5466
                   icond=ivet(4)
5467

    
5468
                case default
5469
                   err_refcodes=.true.
5470
                   ERR_RefCodes_Mess="Only numbers in "//trim(line)
5471
                   return
5472
             end select
5473

    
5474
             nb=0
5475
             do na=1,FAtom%natoms
5476
                call fill_refcodes(key,dire,na,nb,x_low,x_up,x_step,icond,fatom,spg)
5477
             end do
5478
             if (err_refcodes) return
5479

    
5480
          else
5481
             !---- [LABEL, [INF,[SUP,[STEP,[COND]]]]] ----!
5482
             ! If Ilabel(i) == 0 then is a label otherwise is a number
5483
             n_ini=minloc(ilabel,dim=1)
5484
             ilabel(n_ini)=2
5485
             n_end=minloc(ilabel,dim=1)-1
5486

    
5487
             do n=1,ndir
5488
                na=0
5489
                nb=0
5490

    
5491
                !---- Default values ----!
5492
                x_low =0.0
5493
                x_up  =0.0
5494
                x_step=0.0
5495
                icond =0
5496

    
5497
                !---- Label ----!
5498
                npos=index(label(n_ini),"_")
5499
                if (npos >0) then
5500
                   do j=1,ncode
5501
                      if (u_case(label(n_ini)(1:npos))==u_case(trim(code_nam(j)))) then
5502
                         nb=j
5503
                         exit
5504
                      end if
5505
                   end do
5506
                   if (nb == 0) then
5507
                      err_refcodes=.true.
5508
                      ERR_RefCodes_Mess="Code-name not found for "//trim(label(n_ini))
5509
                      return
5510
                   end if