	PROGRAM mk_model

C       F. AFFOUARD 26/11/2014
C	include improper
C	correc bbox 


	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        PARAMETER(N=999999)
        PARAMETER(N_FF=2)
        PARAMETER(NATOM_MAX=100)
        PARAMETER(NSPEC_MAX=9)
        PARAMETER(NOP_MAX=108)
        PARAMETER(CALO=4.187)
        PARAMETER(FAC2=2.0)

        CHARACTER*60 HEADER(N)
        CHARACTER*130 longchain
        CHARACTER*60 COMPOUND_NAME
        CHARACTER*11 MOLECULE_NAME(NSPEC_MAX)
        CHARACTER*13 POTENTIEL_NAME
        DIMENSION NATMS(NATOM_MAX)
        INTEGER SITE_SPEC(NOP_MAX)
        INTEGER TYPE(NATOM_MAX,NSPEC_MAX)
        DIMENSION BOX(3,3)
        DIMENSION INDEX_NB(2,N),POT_NB(3,N)
        DIMENSION TSYM(NOP_MAX,3,4)
        DIMENSION CHARGE(NATOM_MAX,NSPEC_MAX)
        REAL MASS(NATOM_MAX,NSPEC_MAX)
        REAL*8 MASS_ATOM
        DIMENSION X(NATOM_MAX,NSPEC_MAX)
	DIMENSION Y(NATOM_MAX,NSPEC_MAX)
	DIMENSION Z(NATOM_MAX,NSPEC_MAX)
        CHARACTER*7 CHEMICAL(NATOM_MAX,NSPEC_MAX)
        CHARACTER*4 string4
        CHARACTER*8 string8
        CHARACTER*8 GENERATE
	CHARACTER*4 TO_COMPLETE
	CHARACTER*5 FF,READ_FF
	CHARACTER*6 INTER_14
	CHARACTER*40 F_IN,OUT_DLPOLY_CONFIG,OUT_DLPOLY_FIELD
	CHARACTER*40 OUT_SYS_PDB
	CHARACTER*4 READ_ATOM
	CHARACTER*5  INTRA_PROP
	CHARACTER*6  CHAINE
	DIMENSION LISTE_ATOM_RIGID(NATOM_MAX,NSPEC_MAX)
	DIMENSION PARAM(N)
	integer iposchain(4)
        CHARACTER*20 TYPE_TORSION(N)

        CHARACTER*6 CHEM_MMTK(NATOM_MAX,NSPEC_MAX),DUM

	DIMENSION INDEX_CONST(2,N,NSPEC_MAX),POT_CONST(1,N,NSPEC_MAX)
	DIMENSION INDEX_ANGLE(3,N,NSPEC_MAX),POT_ANGLE(2,N,NSPEC_MAX)
	DIMENSION INDEX_BONDS(2,N,NSPEC_MAX),POT_BONDS(2,N,NSPEC_MAX)
        DIMENSION INDEX_TORSION(4,N,NSPEC_MAX)
        DIMENSION INDEX_TORSION_TMP(4,N,NSPEC_MAX)
        DIMENSION POT_TORSION(5,N,NSPEC_MAX)
        DIMENSION INDEX_HARMD(4,N,NSPEC_MAX),POT_HARMD(2,N,NSPEC_MAX)
	DIMENSION INDEX_IMPRO(4,N,NSPEC_MAX),POT_IMPRO(2,N,NSPEC_MAX)

        DIMENSION  N_RIGID(NSPEC_MAX),N_CONST(NSPEC_MAX)
        DIMENSION  N_ANGLE(NSPEC_MAX),N_TORSION(NSPEC_MAX)
	DIMENSION  N_HARMD(NSPEC_MAX),N_BONDS(NSPEC_MAX)
	DIMENSION  N_IMPRO(NSPEC_MAX)

	CHARACTER*12 UNITE
	CHARACTER*11 MOLECULE
	CHARACTER*3 CHEM_TYPE(N)
	CHARACTER*7 DLPOLY_TYPE(N)
        INTEGER VDW_TYPE(N)
	CHARACTER*3 CHEM,CHEM1,CHEM2,CHEM3,CHEM4
	DIMENSION BOX_DIL(3,3)
        INTEGER IVDW(NATOM_MAX,NSPEC_MAX)

	DIMENSION PARA1(N),PARA2(N)
        CHARACTER*5 ORTHO
	CHARACTER*4 DUMMY
	INTEGER TK
	REAL MK
	LOGICAL OK,bigbox
        DIMENSION BOX_INV(3,3)
		
	PI = ACOS(-1.)
	RD = PI/180.
	IDUM = 0


	PRINT *,'make CONFIG & FIELD for DL-POLY'
	PRINT *,'VERSION 2'

	OUT_DLPOLY_CONFIG='CONFIG'
	  PRINT *,'Output dl-poly config : ',OUT_DLPOLY_CONFIG
	OUT_DLPOLY_FIELD='FIELD'
	  PRINT *,'Output dl-poly field  : ',OUT_DLPOLY_FIELD

	DO I=1,NSPEC_MAX
	  N_RIGID(I) = 0
	  N_BONDS(I) = 0
	  N_ANGLE(I) = 0
	  N_TORSION(I) = 0
	ENDDO
	NH = 0	
	OK = .TRUE.

	DO I=1,N
	  IF(OK) THEN 
	    READ(*,'(a60)',end=1) HEADER(I)
	    IF(HEADER(I) .NE. '##' .OR. I.EQ. 1) THEN 
              PRINT *,HEADER(I)
	      NH = NH + 1
	    ENDIF
	    IF(HEADER(I)(1:2) .EQ. '##' .AND. I.GT.1) THEN 
              OK = .FALSE.
	    ENDIF
	  ENDIF
	ENDDO

1       READ(*,'(a,i7)') COMPOUND_NAME
	PRINT *,COMPOUND_NAME

        READ(*,*) A,B,C,ALPHA,BETA,GAMMA,MX,MY,MZ,NGENOP,NSPEC

	  bigbox = .false.
	IF( (MX .EQ. 1) .AND. (MY .EQ. 1) .AND. (MZ .EQ. 1) 
     $      .AND. (NGENOP .EQ. 1) .AND. (NSPEC .EQ. 1)  ) THEN
	  print *,'generate big box 50 A for DLPOLY CONFIG file'
	  bigbox = .true.
	ENDIF

        PRINT '(6a8,5a3)','a','b','c','alpha','beta','gamma',
     $      'mx','my','mz','Z ','spec'
        PRINT '(6F8.3,5I3)',A,B,C,ALPHA,BETA,GAMMA,MX,MY,MZ,NGENOP,NSPEC
       IF(NSPEC .GT. NSPEC_MAX) THEN
          PRINT *,'Augmenter NSPEC_MAX !!!'
          PRINT *,'ne pas oublier les subroutines !!!'
          STOP
        ENDIF

        NMOL = MX*MY*MZ

C       MAKE SIMULATION BOX

        CALL MK_BOX(BOX,A,B,C,ALPHA,BETA,GAMMA)
        CALL INVMAT(BOX,BOX_INV)

        READ(*,'(a)') ORTHO
        IF(ORTHO .NE. 'fract' .AND. ORTHO .NE. 'ortho') THEN
          PRINT *,'ORTHO = fract or ortho !!!'
          STOP
        ENDIF

	PRINT *,'Box = '
	WRITE(*,'(3F10.5)') BOX(1,1),BOX(1,2),BOX(1,3)
	WRITE(*,'(3F10.5)') BOX(2,1),BOX(2,2),BOX(2,3)
	WRITE(*,'(3F10.5)') BOX(3,1),BOX(3,2),BOX(3,3)
	PRINT *,' '

	
C       FORCE_FIELD
        READ(*,'(a)') FF
        OPEN(UNIT=14,FILE='/mnt/c/bin/FF/ff.list',STATUS='OLD')
        CHK_FF = 0
        DO IFF=1,N_FF
          READ(14,'(a)') READ_FF
          IF(READ_FF .EQ. FF) THEN
            PRINT *,'Force Field = ',FF
            CHK_FF = 1
          ENDIF
        ENDDO
        IF(CHK_FF .EQ. 0.0) THEN
          PRINT *,'!!! unknown force field : ',FF
          PRINT *,'increase N_FF or check force field'
          STOP
        ENDIF
        CLOSE(14)

	READ(*,*) INTER_14,RESCALE_VDW_14,RESCALE_ELEC_14

	IF(INTER_14(1:3) .EQ. 'e14') THEN
	  PRINT '(a,2F5.2)','Rescale 1-4 interaction factor = ',
     $  RESCALE_VDW_14,RESCALE_ELEC_14
	ELSE
	  PRINT *,'rescale 1-4 factor : simple '
	  STOP
	ENDIF

        UNITE = 'kcal'
        PRINT '(2a)','... Physical Unit Read = ',UNITE
        CAL2J = 1.0

	PRINT *,' '
        ii = 0
        CHEM = ' '

C	LIRE ATOMIC SPECIFICATION
	IK = 0
        DO J=1,NSPEC
	  READ(*,*) MOLECULE_NAME(J),NATMS(J)

C	DEBUG
        OK = .TRUE.
	L = 11
        DO IL=1,11
          IF(OK .AND. MOLECULE_NAME(J)(IL:IL) .EQ. ' ') THEN
            L = IL - 1
            OK = .FALSE.
          ENDIF
        ENDDO
 
	  PRINT *,'SPEC : ',MOLECULE_NAME(J),' ',NATMS(J)
	  PRINT '(a5,5a9)','atom','x','y','z','mass','q'

	  CHARGETOT = 0.
          DO IA=1,NATMS(J)

            READ(*,*)
     $      CHEM_MMTK(IA,J), III,XK,YK,ZK,READ_Q,READ_ATOM,IREAD_VDW

	    CHARGE(IA,J)   = READ_Q
            IVDW(IA,J)  = IREAD_VDW
            
            I4 = IREAD_VDW/1000
            I3 = (IREAD_VDW - I4*1000)/100
            I2 = (IREAD_VDW - I4*1000 - I3*100)/10
            I1 = IREAD_VDW - I4*1000 - I3*100 -I2*10
            string4(1:1) = CHAR(48+I4)
            string4(2:2) = CHAR(48+I3)
            string4(3:3) = CHAR(48+I2)
            string4(4:4) = CHAR(48+I1)
            string8 = READ_ATOM//string4
            call cleanstring(string8)
            CHEMICAL(IA,J) = string8

            CALL 
     $      FF_PARAM(FF,'masse',READ_ATOM,CHEM,CHEM,CHEM,
     $      1,IREAD_VDW,PARAM,longchain,TYPE_TORSION)
            READ_MASS      = PARAM(1)
            MASS(IA,J)     = READ_MASS
          
C       cherche position "
          k  = 0
          do i=1,130
            if(longchain(i:i) .eq. '"') then 
              k = k + 1
              if(k .gt. 4) then
                 print *,'!! problem file masse.tab'
                 print *,'too many " '
                 stop
              endif 
              iposchain(k) = i
            endif 
          enddo
 

            CALL 
     $      FF_PARAM(FF,'vdw1a',READ_ATOM,CHEM,CHEM,CHEM,
     $      3,IREAD_VDW,PARAM,longchain,TYPE_TORSION)
            QQ   =  PARAM(1)


            if(QQ .ne. READ_Q) then
               PRINT *,'!!!!'
               PRINT '(2i5,1x,a8)',IA,IVDW(IA,J),CHEMICAL(IA,J) 
               print *,'charges do not match'
               print *,'in .def file = ',READ_Q
               print *,'in .tab file = ',QQ
            endif

           IF(IREAD_VDW .NE. ii) THEN
              OK = .TRUE.
              INEW = 0.0
              DO ITEST=1,IK
                IF(IREAD_VDW .EQ. VDW_TYPE(ITEST)) THEN
                  OK = .FALSE.
                  IGET = ITEST
                ENDIF
              ENDDO ! ITEST
              TYPE(IA,J) = INEW
              IF( OK ) THEN
                IK = IK + 1
                VDW_TYPE(IK) = IREAD_VDW
                CHEM_TYPE(IK) = READ_ATOM
                DLPOLY_TYPE(IK) = CHEMICAL(IA,J)
                IGET = IK
              ENDIF
            ENDIF
            TYPE(IA,J) = IGET
            ii = IREAD_VDW

	    IF(ORTHO .EQ. 'ortho') 
     $      CALL ORTH2CRYS(XK,YK,ZK,BOX_INV)

	    X(IA,J)    = XK
	    Y(IA,J)    = YK 
	    Z(IA,J)    = ZK


C	CONVERSION CARTESIAN COORD
	    IF(ORTHO .NE. 'ortho') 
     $      CALL ORTH2CRYS(XK,YK,ZK,BOX)

            PRINT '(2i4,2a5,5F9.4,a,1x,a,1x,a,1x,a)',IA,IVDW(IA,J),
     $          READ_ATOM,CHEM_MMTK(IA,J),XK,YK,ZK, 
     $          READ_MASS,READ_Q,' ',CHEMICAL(IA,J),
     $              longchain(iposchain(1):iposchain(2)),
     $              longchain(iposchain(3):iposchain(4))       
	
	    CHARGETOT = CHARGETOT + READ_Q
          ENDDO ! NATMS
	  IF(ABS(CHARGETOT) .GT. 0.00001) THEN
	    print '(a20,f10.5)', 'Charge totale = ',CHARGETOT
	    STOP
	  ENDIF
	  PRINT *,' '
          ENDDO ! NSPECIES


          N_TYPE  = IK
          PRINT *,'N different chemicals = ',N_TYPE
          PRINT *,' '

	  PRINT *,' '
	  PRINT '(a12,3a8)','Species','x','y','z'
	  PRINT '(3a8)','SymOp_x','SymOp_y','SymOp_z'

	NOP = NGENOP*NSPEC
        IF(NOP .GT. NOP_MAX) THEN
          PRINT *,'Augmenter NOP_MAX !!! '
          PRINT *,'ne pas oublier les subroutines !!!'
          STOP
        ENDIF

	DO J2=1,NSPEC
          DO J1=1,NGENOP
	    J = J1 + NGENOP*(J2-1)
	    SITE_SPEC(J) = J2
	    IF(J2 .EQ. 1) THEN 
              READ(*,*) TSYM(J,1,4),TSYM(J,2,4),TSYM(J,3,4)
              READ(*,*) TSYM(J,1,1),TSYM(J,1,2),TSYM(J,1,3)
              READ(*,*) TSYM(J,2,1),TSYM(J,2,2),TSYM(J,2,3)
              READ(*,*) TSYM(J,3,1),TSYM(J,3,2),TSYM(J,3,3)
              PRINT '(i5,a12,3F8.2)',J1,MOLECULE_NAME(SITE_SPEC(J)),
     $        TSYM(J,1,4),TSYM(J,2,4),TSYM(J,3,4)
              PRINT '(3F8.2)', TSYM(J,1,1),TSYM(J,1,2),TSYM(J,1,3)
              PRINT '(3F8.2)', TSYM(J,2,1),TSYM(J,2,2),TSYM(J,2,3)
              PRINT '(3F8.2)', TSYM(J,3,1),TSYM(J,3,2),TSYM(J,3,3)
	    ELSE
	      DO K1=1,3
              DO K2=1,4
                TSYM(J,K1,K2) = TSYM(J1,K1,K2)
	      ENDDO
	      ENDDO
	    ENDIF
          ENDDO !NGENOP
        ENDDO !NSPEC

	PRINT *,' '

C	DIPOLE MOMENT 
	
	PRINT *,'Dipole moment (Debye) : '
        CALL CAL_DIPOLE(NATOM_MAX,NSPEC_MAX,NSPEC,
     $  NATMS,X,Y,Z,BOX,CHARGE,MOLECULE_NAME)

C	INTRAMOLECULAR INTERACTION

	DO ISPEC=1,NSPEC

	  READ(*,*) MOLECULE,NINTRA
	  PRINT *,' '
	  PRINT *,'*** Intramolecular interaction for spec = ',MOLECULE
          PRINT *,NINTRA,' interaction(s) expected'
	  
          IF( MOLECULE_NAME(ISPEC) .NE. MOLECULE) THEN
             PRINT *,'spec expected = ',MOLECULE_NAME(ISPEC),'!!!'
             STOP
          ENDIF
	
	  DO INTRA=1,NINTRA

	  READ(*,*) INTRA_PROP,GENERATE,NP

	  IF(GENERATE .EQ. 'generate') THEN 
            PRINT '(a5,1x,a8)',INTRA_PROP,GENERATE
	  ELSE
            PRINT '(a5,1x,I5)',INTRA_PROP,NP
	  ENDIF

C RIGID
	  IF(INTRA_PROP .EQ. 'rigid') THEN
	    N_RIGID(ISPEC) = NP
	    DO J=1,NP
	      LISTE_ATOM_RIGID(J,ISPEC) = J
	    ENDDO
	  ENDIF

C CONST
          IF(INTRA_PROP .EQ. 'const') THEN
            N_CONST(ISPEC)=NP
            PRINT *,'constraint bond'
            PRINT *,'(1/2 for DL_POLY !!!)'
            PRINT '(2a5,2a5,a12)','atom','atom','chem','chem','r0'
            DO J=1,NP
              READ(*,*) (INDEX_CONST(K,J,ISPEC),K=1,2),
     $                POT_CONST(1,J,ISPEC)
              I1 = INDEX_CONST(1,J,ISPEC)
              I2 = INDEX_CONST(2,J,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              PRINT '(2i5,2a5,f12.5)',INDEX_CONST(1,J,ISPEC),
     $                   INDEX_CONST(2,J,ISPEC),CHEM1,CHEM2,
     $                   POT_CONST(1,J,ISPEC) 
            ENDDO
          ENDIF


C BONDS
          IF(INTRA_PROP .EQ. 'bonds') THEN
            N_BONDS(ISPEC)=NP
	    PRINT *,'V(r) = 1/2k(r-r0)^2'
	    PRINT *,'(1/2 for DL_POLY !!!)'
	    PRINT '(2a5,2a5,2a12)','atom','atom','chem','chem','k','r0'
            DO J=1,NP
              READ(*,*) (INDEX_BONDS(K,J,ISPEC),K=1,2)
	      I1 = INDEX_BONDS(1,J,ISPEC)
	      I2 = INDEX_BONDS(2,J,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)

              CALL 
     $        FF_PARAM(FF,INTRA_PROP,CHEM1,CHEM2,CHEM,CHEM,
     $        1,IDUM,PARAM,longchain,TYPE_TORSION)
              POT_BONDS(1,J,ISPEC) = CAL2J*PARAM(1)
              POT_BONDS(2,J,ISPEC) = PARAM(2)


	      PRINT '(2i5,2a5,2f12.5)',INDEX_BONDS(1,J,ISPEC),
     $                   INDEX_BONDS(2,J,ISPEC),CHEM1,CHEM2,
     $                   POT_BONDS(1,J,ISPEC),
     $                   POT_BONDS(2,J,ISPEC)
            ENDDO
          ENDIF



C ANGLE

	  IF(INTRA_PROP .EQ. 'angle') THEN
            PRINT *,'V(th) = 1/2k(th-th0)^2'
	    PRINT *,'(1/2 for DL_POLY !!!)'
            PRINT '(3a5,3a5,2a12)','atom','atom','atom',
     $                             'chem','chem','chem','k','th0'

	    IF(GENERATE .EQ. 'generate') THEN 
              PRINT *,'...generate angles list fron bonds list'
              IF(N_BONDS(ISPEC) .EQ. 0 ) THEN
                PRINT *,'bonds list or constraint list are required'
                STOP
              ENDIF
                CALL GENERATE_ANGLE_LISTE(ISPEC,N,NSPEC_MAX,
     $          INDEX_BONDS,N_BONDS(ISPEC),NP,INDEX_ANGLE,NATMS(ISPEC))
	    ELSE
	      DO J=1,NP
                READ(*,*) (INDEX_ANGLE(K,J,ISPEC),K=1,3)
	      ENDDO
	    ENDIF
	    N_ANGLE(ISPEC) = NP

	    DO J=1,NP
              I1 = INDEX_ANGLE(1,J,ISPEC)
              I2 = INDEX_ANGLE(2,J,ISPEC)
              I3 = INDEX_ANGLE(3,J,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              IK3 = TYPE(I3,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              CHEM3 = CHEM_TYPE(IK3)
              X1 =  X(I1,ISPEC)
              Y1 =  Y(I1,ISPEC)
              Z1 =  Z(I1,ISPEC)
              X2 =  X(I2,ISPEC)
              Y2 =  Y(I2,ISPEC)
              Z2 =  Z(I2,ISPEC)
              X3 =  X(I3,ISPEC)
              Y3 =  Y(I3,ISPEC)
              Z3 =  Z(I3,ISPEC)
              CALL
     $        FF_PARAM(FF,INTRA_PROP,CHEM1,CHEM2,CHEM3,CHEM,
     $        1,IDUM,PARAM,longchain,TYPE_TORSION)
              POT_ANGLE(1,J,ISPEC) = PARAM(1)
              POT_ANGLE(2,J,ISPEC) = PARAM(2)

              PRINT '(3i5,3a5,2f12.5)',
     $    INDEX_ANGLE(1,J,ISPEC),
     $    INDEX_ANGLE(2,J,ISPEC),
     $    INDEX_ANGLE(3,J,ISPEC),
     $    CHEM1,CHEM2,CHEM3,
     $    POT_ANGLE(1,J,ISPEC),POT_ANGLE(2,J,ISPEC)
	    ENDDO
	  ENDIF


C IMPRO
	IF(INTRA_PROP .EQ. 'impro') THEN
           PRINT *,'...read improper list '
           DO J=1,NP
             READ(*,*) (INDEX_IMPRO(K,J,ISPEC),K=1,4)
            ENDDO
           N_IMPRO(ISPEC) = NP
        ENDIF

C TORSI

          IF(INTRA_PROP .EQ. 'torsi') THEN

C       read intra type torsion

             read(*,*) longchain,n_type_torsion
             print '(a,i5)','n type torsion = ',n_type_torsion
             do kk=1,n_type_torsion
               read(*,'(a20)') type_torsion(kk)
               print *,kk,' torsion ok = ',type_torsion(kk)
             enddo 

 	    PRINT *,'V(r) = k[1+cos(m.th-delta)]'
            PRINT '(4a5,4a5,3a12)','atom','atom','atom','atom',
     $                             'chem','chem','chem','chem',
     $                             'k','m','delta'

            IF(GENERATE .EQ. 'generate') THEN
              PRINT *,'...generate torsion list from bonds list'
              IF(N_BONDS(ISPEC) .EQ. 0 .AND.
     $          N_CONST(ISPEC) .EQ. 0) THEN
                PRINT *,'bonds list or constraint list is required'
                STOP
              ENDIF
              IF(N_CONST(ISPEC) .EQ. 0) THEN
                CALL GENERATE_TORSI_LISTE(ISPEC,N,NSPEC_MAX,
     $          INDEX_BONDS,N_BONDS(ISPEC),NP,INDEX_TORSION_TMP
     $          ,NATMS(ISPEC))
              ELSE
                CALL GENERATE_TORSI_LISTE(ISPEC,N,NSPEC_MAX,
     $          INDEX_CONST,N_CONST(ISPEC),NP,INDEX_TORSION_TMP
     $          ,NATMS(ISPEC))
              ENDIF
            ELSE
              PRINT *,'...read torsion list '
              DO J=1,NP
                  READ(*,*) (INDEX_TORSION_TMP(K,J,ISPEC),K=1,4)
              ENDDO
            ENDIF
            N_TORSION(ISPEC) = NP

	    J = 0
            DO JJ=1,NP
              I1 = INDEX_TORSION_TMP(1,JJ,ISPEC)
              I2 = INDEX_TORSION_TMP(2,JJ,ISPEC)
              I3 = INDEX_TORSION_TMP(3,JJ,ISPEC)
              I4 = INDEX_TORSION_TMP(4,JJ,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              IK3 = TYPE(I3,ISPEC)
              IK4 = TYPE(I4,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              CHEM3 = CHEM_TYPE(IK3)
              CHEM4 = CHEM_TYPE(IK4)
              X1 =  X(I1,ISPEC)
              Y1 =  Y(I1,ISPEC)
              Z1 =  Z(I1,ISPEC)
              X2 =  X(I2,ISPEC)
              Y2 =  Y(I2,ISPEC)
              Z2 =  Z(I2,ISPEC)
              X3 =  X(I3,ISPEC)
              Y3 =  Y(I3,ISPEC)
              Z3 =  Z(I3,ISPEC)
              X4 =  X(I4,ISPEC)
              Y4 =  Y(I4,ISPEC)
              Z4 =  Z(I4,ISPEC)


              P4 = RESCALE_VDW_14
	      P5 = RESCALE_ELEC_14


	      IF(FF .EQ. 'OPLS_') THEN 

C               V1,V2,V3
C               P1,P2,P3
                  P1 = 0.
	          P2 = 0.
                  P3 = 0. 

                  CALL
     $            FF_PARAM(FF,INTRA_PROP,CHEM1,CHEM2,
     $            CHEM3,CHEM4,Imatch,n_type_torsion,PARAM,
     $            longchain,TYPE_TORSION)

C       only Imatch = 1 is ok

                  IF(Imatch .eq. 1) THEN 

                  P1 = CAL2J*PARAM(1)
                  P2 = PARAM(2)
                  P3 = PARAM(3)
	          P1 = 0.5*P1
	          P2 = 0.5*P2
	          P3 = 0.5*P3
                  IF(P1 .NE. 0.) THEN
                    J = J + 1
                    INDEX_TORSION(1,J,ISPEC) = I1
                    INDEX_TORSION(2,J,ISPEC) = I2
                    INDEX_TORSION(3,J,ISPEC) = I3
                    INDEX_TORSION(4,J,ISPEC) = I4
                    POT_TORSION(1,J,ISPEC) = P1
                    POT_TORSION(2,J,ISPEC) = 1.0
                    POT_TORSION(3,J,ISPEC) = 0.0 
                    POT_TORSION(4,J,ISPEC) = P4
                    POT_TORSION(5,J,ISPEC) = P5
                  ENDIF
                  IF(P2 .NE. 0.) THEN
                    J = J + 1
                    INDEX_TORSION(1,J,ISPEC) = I1
                    INDEX_TORSION(2,J,ISPEC) = I2
                    INDEX_TORSION(3,J,ISPEC) = I3
                    INDEX_TORSION(4,J,ISPEC) = I4
                    POT_TORSION(1,J,ISPEC) = P2
                    POT_TORSION(2,J,ISPEC) = 2.0
                    POT_TORSION(3,J,ISPEC) = 180.0
                    POT_TORSION(4,J,ISPEC) = P4
                    POT_TORSION(5,J,ISPEC) = P5
                  ENDIF
                  IF(P3 .NE. 0.) THEN
                    J = J + 1
                    INDEX_TORSION(1,J,ISPEC) = I1
                    INDEX_TORSION(2,J,ISPEC) = I2
                    INDEX_TORSION(3,J,ISPEC) = I3
                    INDEX_TORSION(4,J,ISPEC) = I4
                    POT_TORSION(1,J,ISPEC) = P3
                    POT_TORSION(2,J,ISPEC) = 3.0
                    POT_TORSION(3,J,ISPEC) = 0.0
                    POT_TORSION(4,J,ISPEC) = P4
                    POT_TORSION(5,J,ISPEC) = P5
                  ENDIF
                
                endif ! Imatch = 1 

                N_TORSION(ISPEC) = J

	      ENDIF ! 

            ENDDO ! JJ

            N_TORSION(ISPEC) = J
        
            DO J=1,N_TORSION(ISPEC)

              I1=INDEX_TORSION(1,J,ISPEC) 
              I2=INDEX_TORSION(2,J,ISPEC) 
              I3=INDEX_TORSION(3,J,ISPEC)
              I4=INDEX_TORSION(4,J,ISPEC) 
              P1=POT_TORSION(1,J,ISPEC) 
              P2=POT_TORSION(2,J,ISPEC)
              P3=POT_TORSION(3,J,ISPEC) 
              P4=POT_TORSION(4,J,ISPEC)
              P5=POT_TORSION(5,J,ISPEC) 

              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              IK3 = TYPE(I3,ISPEC)
              IK4 = TYPE(I4,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              CHEM3 = CHEM_TYPE(IK3)
              CHEM4 = CHEM_TYPE(IK4)

              PRINT '(4i5,4a5,5f12.5)',
     $        I1,I2,I3,I4,
     $        CHEM1,CHEM2,CHEM3,CHEM4,
     $        P1,P2,P3,P4,P5

            ENDDO 
	  print *," total torsions found    = ",NP
	  print *," torsions removed (zero) = ",NP-J


          ENDIF ! opls 


	ENDDO ! INTRAMOLECULAR PROPERTIES
	ENDDO ! NSPEC


C NON BONDED POTENTIAL 

	IF (FF .EQ. 'GAFF_') THEN

          PRINT *,'GAFF (AMBER) non-bonded potential '
          PRINT *,'convert to : V(r) = 4eps.[(sig/r)^12 -(sig/r)^6]'
          PRINT '(2a5,2a12)',
     $ 'chem','chem','eps(LJ)','sig(LJ)'

          print *,'to be impleted for GAFF'
        stop 

        ELSE IF(FF .EQ. 'OPLS_') THEN
          PRINT *,'OPLS non-bonded potential '
          PRINT *,'convert to : V(r) = 4eps.[(sig/r)^12 -(sig/r)^6]'
          PRINT '(2a5,2a12)',
     $ 'chem','chem','eps(LJ)','sig(LJ)'
          J = 0
          DO I1=1,N_TYPE
            DO I2=I1,N_TYPE
              IVDW1 = VDW_TYPE(I1)
              IVDW2 = VDW_TYPE(I2)
              CALL
     $        FF_PARAM(FF,'vdw1a',CHEM1,CHEM,CHEM,CHEM,
     $        1,IVDW1,PARAM,longchain,TYPE_TORSION)
              P11 = PARAM(2)
              P12 = PARAM(3)
              CALL
     $        FF_PARAM(FF,'vdw1a',CHEM2,CHEM,CHEM,CHEM,
     $        1,IVDW2,PARAM,longchain,TYPE_TORSION)
              P21 = PARAM(2)
              P22 = PARAM(3)

              P1  = SQRT(P11*P21)
              P2  = SQRT(P12*P22)
              EPS  = P2
              SIG  = P1

              IF(P1 .NE. 0. .AND. P2 .NE. 0.) THEN
                J = J + 1
                INDEX_NB(1,J) = I1
                INDEX_NB(2,J) = I2
                POT_NB(1,J) = P2
                POT_NB(2,J) = P1
                PRINT '(2a7,2F12.5)',DLPOLY_TYPE(I1),DLPOLY_TYPE(I2),
     $          POT_NB(2,J),POT_NB(1,J)
              ENDIF
            ENDDO
          ENDDO

        
          N_INTERACTION = J
          PRINT *,'N interactions = ',J

	ELSE
          PRINT *,'unknown non-bonded force field !!!'
          STOP
        ENDIF
	

	  OK = .FALSE.
	  DO ISPEC=1,NSPEC
	    IF( (N_CONST(ISPEC) .GT. 0) .OR.
     $          (N_BONDS(ISPEC) .GT. 0) .OR.
     $          (N_ANGLE(ISPEC) .GT. 0) .OR.
     $          (N_HARMD(ISPEC) .GT. 0) .OR.
     $          (N_TORSION(ISPEC) .GT. 0)) OK = .TRUE.
	  ENDDO

	  IF(OK) CALL  CHECK_GEOM(X,Y,Z,NSPEC,
     $  N_CONST,N_BONDS,N_ANGLE,N_HARMD,N_TORSION,
     $                       N_IMPRO,
     $                       INDEX_CONST,POT_CONST,
     $                       INDEX_ANGLE,POT_ANGLE,
     $                       INDEX_BONDS,POT_BONDS,
     $                       INDEX_TORSION,POT_TORSION,
     $                       INDEX_HARMD,POT_HARMD,
     $                       INDEX_IMPRO,POT_IMPRO,
     $                       BOX,NATMS,TYPE,CHEM_TYPE)


	CALL MK_DLPOLY(MX,MY,MZ,NSPEC,NOP,NMOL,N_INTERACTION,
     $  OUT_DLPOLY_CONFIG,OUT_SYS_PDB,OUT_DLPOLY_FIELD,
     $  COMPOUND_NAME,MOLECULE_NAME,POTENTIEL_NAME,
     $  BOX,BOX_DIL,DILATATION,X,Y,Z,TYPE,NATMS,
     $  CHEMICAL,SITE_SPEC,
     $  MASS,CHARGE,INDEX_NB,POT_NB,TSYM,CHEM_TYPE,N_RIGID,
     $  N_CONST,N_BONDS,N_ANGLE,N_HARMD,N_TORSION,N_IMPRO,
     $                       LISTE_ATOM_RIGID,
     $                       INDEX_CONST,POT_CONST,
     $                       INDEX_ANGLE,POT_ANGLE,
     $                       INDEX_BONDS,POT_BONDS,
     $                       INDEX_TORSION,POT_TORSION,
     $                       INDEX_IMPRO,POT_IMPRO,
     $                       INDEX_HARMD,POT_HARMD,FF,
     $                       DLPOLY_TYPE,bigbox)



	PRINT *,'...termine'
	STOP
	END


C**************************************************

	SUBROUTINE MK_BOX(BOX,A,B,C,ALPHA,BETA,GAMMA)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	DIMENSION BOX(3,3)

C       CHANGEMENT REP Xtal coord -> ortho coord
C       A // x & B in (xoy) plane

	PI = DACOS(-1.0D0)
	RD = PI/180.D0
	
	CGAMMA = DCOS(RD*GAMMA)
	SGAMMA = DSIN(RD*GAMMA)
	CBETA  = DCOS(RD*BETA)
	CALPHA = DCOS(RD*ALPHA)

C	BOX(1,1) = A
C	BOX(1,2) = 0.0
C	BOX(1,3) = 0.0
C	BOX(2,1) = B*CGAMMA
C	BOX(2,2) = B*SGAMMA
C	BOX(2,3) = 0.0
C	BOX(3,1) = C*CBETA
C	BOX(3,2) = C*(CALPHA - CGAMMA*CBETA)/SGAMMA
C	V =C*SQRT(1.- CGAMMA**2 -CBETA**2 -CALPHA**2 +2.*CGAMMA*CALPHA*CBETA)
c       BOX(3,3) = V/(SGAMMA)

 	BOX(1,1) = A
 	BOX(2,1) = 0.0
 	BOX(3,1) = 0.0
        BOX(1,2) = B*CGAMMA
        BOX(2,2) = B*SGAMMA
        BOX(3,2) = 0.0
        BOX(1,3) = C*CBETA
        BOX(2,3) = C*(CALPHA - CGAMMA*CBETA)/SGAMMA
 	V =C*SQRT(1.- CGAMMA**2 -CBETA**2 -
     $  CALPHA**2 +2.*CGAMMA*CALPHA*CBETA)
        BOX(3,3) = V/(SGAMMA)
 
	RETURN
	END


C**************************************************

C       CHECK GEOM

	SUBROUTINE CHECK_GEOM(X,Y,Z,NSPEC,
     $  N_CONST,N_BONDS,N_ANGLE,N_HARMD,N_TORSION,
     $                       N_IMPRO,
     $                       INDEX_CONST,POT_CONST,
     $                       INDEX_ANGLE,POT_ANGLE,
     $                       INDEX_BONDS,POT_BONDS,
     $                       INDEX_TORSION,POT_TORSION,
     $                       INDEX_HARMD,POT_HARMD,
     $                       INDEX_IMPRO,POT_IMPRO,
     $                       BOX,NATMS,TYPE,CHEM_TYPE)

        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        PARAMETER(NATOM_MAX=100)
        PARAMETER(NSPEC_MAX=9)
        PARAMETER(N=999999)
        PARAMETER(EPSILON=.10)

        DIMENSION X(NATOM_MAX,NSPEC_MAX)
        DIMENSION Y(NATOM_MAX,NSPEC_MAX)
        DIMENSION Z(NATOM_MAX,NSPEC_MAX)
	DIMENSION XA(NATOM_MAX)
	DIMENSION YA(NATOM_MAX)
	DIMENSION ZA(NATOM_MAX)
	DIMENSION BOX(3,3),NATMS(NATOM_MAX)

        DIMENSION INDEX_CONST(2,N,NSPEC_MAX),POT_CONST(1,N,NSPEC_MAX)
        DIMENSION INDEX_ANGLE(3,N,NSPEC_MAX),POT_ANGLE(2,N,NSPEC_MAX)
        DIMENSION INDEX_BONDS(2,N,NSPEC_MAX),POT_BONDS(2,N,NSPEC_MAX)
        DIMENSION INDEX_TORSION(4,N,NSPEC_MAX)
        DIMENSION POT_TORSION(5,N,NSPEC_MAX)
        DIMENSION INDEX_HARMD(4,N,NSPEC_MAX),POT_HARMD(2,N,NSPEC_MAX)
	DIMENSION INDEX_IMPRO(4,N,NSPEC_MAX),POT_IMPRO(2,N,NSPEC_MAX)

	CHARACTER*3 CHEM_TYPE(N),CHEM1,CHEM2,CHEM3,CHEM4
	INTEGER TYPE(NATOM_MAX,NSPEC_MAX)

        DIMENSION  N_RIGID(NSPEC_MAX),N_CONST(NSPEC_MAX)
        DIMENSION  N_ANGLE(NSPEC_MAX),N_TORSION(NSPEC_MAX)
        DIMENSION  N_HARMD(NSPEC_MAX),N_BONDS(NSPEC_MAX)
	DIMENSION  N_IMPRO(NSPEC_MAX)

	PI = ACOS(-1.)
	RD = 180./PI

	PRINT *,'check molecular geometry'

	DO ISPEC=1,NSPEC

          DO K=1,NATMS(ISPEC)
            XA(K) = X(K,ISPEC)
            YA(K) = Y(K,ISPEC)
            ZA(K) = Z(K,ISPEC)
	  ENDDO

          IF(N_CONST(ISPEC) .GT. 0) THEN
	    PRINT *,'...check constraint'
            DO J=1,N_CONST(ISPEC)
              I1 = INDEX_CONST(1,J,ISPEC)
              I2 = INDEX_CONST(2,J,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              X1 =  XA(I1)
              Y1 =  YA(I1)
              Z1 =  ZA(I1)
              X2 =  XA(I2)
              Y2 =  YA(I2)
              Z2 =  ZA(I2)
              R1  = BONDS(X1,Y1,Z1,X2,Y2,Z2,BOX)
	      EXPECT= POT_CONST(1,J,ISPEC)
	      DIST = ABS(R1 - EXPECT)/EXPECT
              IF(DIST .GT. EPSILON) 
     $        PRINT '(2i5,2a5,2f12.5)',I1,I2,CHEM1,CHEM2,
     $                   EXPECT,R1

            ENDDO
          ENDIF


          IF(N_BONDS(ISPEC) .GT. 0) THEN
	    PRINT *,'...check bonds'
              PRINT '(2a5,2a5,2a12)'," "," ",
     $        " "," ","ff","real"
            DO J=1,N_BONDS(ISPEC)
              I1 = INDEX_BONDS(1,J,ISPEC)
              I2 = INDEX_BONDS(2,J,ISPEC)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              X1 =  XA(I1)
              Y1 =  YA(I1)
              Z1 =  ZA(I1)
              X2 =  XA(I2)
              Y2 =  YA(I2)
              Z2 =  ZA(I2)

              R1  = BONDS(X1,Y1,Z1,X2,Y2,Z2,BOX)

              EXPECT= POT_BONDS(2,J,ISPEC)

              IF(EXPECT .EQ. 0.0) THEN
                DIST = ABS(R1 - EXPECT)
              ELSE
                DIST = ABS(R1 - EXPECT)/EXPECT
              ENDIF

              IF(DIST .GT. EPSILON)
     $        PRINT '(2i5,2a5,2f12.5)',I1,I2,CHEM1,CHEM2,
     $                   EXPECT,R1

            ENDDO
          ENDIF

          IF(N_ANGLE(ISPEC) .GT. 0) THEN
	    PRINT *,'...check angles'
              PRINT '(3a5,3a5,2a12)'," "," "," ",
     $        " "," "," ","ff","real"
            DO J=1,N_ANGLE(ISPEC)
              I1 = INDEX_ANGLE(1,J,ISPEC)
              I2 = INDEX_ANGLE(2,J,ISPEC)
              I3 = INDEX_ANGLE(3,J,ISPEC)
C	      CHECK IF BOND(I1-I2) EXIST
C	      CHECK IF BOND(I3-I2) EXIST

	      CALL CHECK_IF_BOND_EXIST(I1,I2,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)

	      CALL CHECK_IF_BOND_EXIST(I3,I2,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)

              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              IK3 = TYPE(I3,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              CHEM3 = CHEM_TYPE(IK3)
	
              X1 =  XA(I1)
              Y1 =  YA(I1)
              Z1 =  ZA(I1)
              X2 =  XA(I2)
              Y2 =  YA(I2)
              Z2 =  ZA(I2)
              X3 =  XA(I3)
              Y3 =  YA(I3)
              Z3 =  ZA(I3)

	      R1 =  ANGLE(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,BOX)

              EXPECT= POT_ANGLE(2,J,ISPEC)
              IF(EXPECT .EQ. 0.0) THEN 
                DIST = ABS(R1 - EXPECT)
              ELSE
                DIST = ABS(R1 - EXPECT)/EXPECT
              ENDIF

	         IF(DIST .GT. EPSILON)
     $        PRINT '(3i5,3a5,2f12.5)',I1,I2,I3,
     $        CHEM1,CHEM2,CHEM3,EXPECT,R1

            ENDDO
          ENDIF


          IF(N_TORSION(ISPEC) .GT. 0) THEN

            PRINT *,'...check torsions'
              PRINT '(4a5,4a5,2a12)'," "," "," "," ",
     $        " "," "," "," ","ff","real"
            DO J=1,N_TORSION(ISPEC)
              I1 = INDEX_TORSION(1,J,ISPEC)
              I2 = INDEX_TORSION(2,J,ISPEC)
              I3 = INDEX_TORSION(3,J,ISPEC)
              I4 = INDEX_TORSION(4,J,ISPEC)
	      CALL CHECK_IF_BOND_EXIST(I1,I2,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)
	      CALL CHECK_IF_BOND_EXIST(I2,I3,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)
	      CALL CHECK_IF_BOND_EXIST(I3,I4,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)
              IK1 = TYPE(I1,ISPEC)
              IK2 = TYPE(I2,ISPEC)
              IK3 = TYPE(I3,ISPEC)
              IK4 = TYPE(I4,ISPEC)
              CHEM1 = CHEM_TYPE(IK1)
              CHEM2 = CHEM_TYPE(IK2)
              CHEM3 = CHEM_TYPE(IK3)
              CHEM4 = CHEM_TYPE(IK4)

              X1 =  XA(I1)
              Y1 =  YA(I1)
              Z1 =  ZA(I1)
              X2 =  XA(I2)
              Y2 =  YA(I2)
              Z2 =  ZA(I2)
              X3 =  XA(I3)
              Y3 =  YA(I3)
              Z3 =  ZA(I3)
              X4 =  XA(I4)
              Y4 =  YA(I4)
              Z4 =  ZA(I4)

	      R1 = DIHED(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,BOX)
              PAR2= POT_TORSION(2,J,ISPEC)
              PAR3= POT_TORSION(3,J,ISPEC)

C	extrema = 
	      EXPECT = PAR2/PAR3 
C       cherche min 
	      KP = 0
	      
	      IF (COS(EXPECT/RD) .GT. COS(EXPECT/RD+PI/PAR3)) THEN
                KP = 1
	      ENDIF
	      EXPECT1 = EXPECT + REAL(KP)*(PI*RD)/PAR3
	      EXPECT2 = EXPECT - REAL(KP)*(PI*RD)/PAR3
	      IF( EXPECT1 .GE. EXPECT2) THEN 
                EXPECT = EXPECT2
	      ELSE
                EXPECT = EXPECT1
	      ENDIF

              IF(EXPECT .EQ. 0.0) THEN
                DIST = ABS(R1 - EXPECT)
              ELSE
                DIST = ABS(R1 - EXPECT)/EXPECT
              ENDIF

              IF(DIST .GT. EPSILON) THEN 
              PRINT '(4i5,4a5,2f12.5)',I1,I2,I3,I4,
     $        CHEM1,CHEM2,CHEM3,CHEM4,EXPECT,R1
	      ENDIF 
 
            ENDDO


          ENDIF

        ENDDO ! NSPEC

	RETURN
	END



C**************************************************

        SUBROUTINE QUATERNION(Q0,Q1,Q2,Q3,TSYM,BOX,IOP)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	PARAMETER(NOP_MAX = 108)
	DIMENSION TSYM(NOP_MAX,3,4)
	DIMENSION EPS(3,3)
	DIMENSION BOX(3,3)
	DIMENSION R(3,3)
	DIMENSION R1(3,3)
	REAL*8      INVBOX(3,3)


	DO I=1,3
	  DO J=1,3
	    EPS(I,J) = TSYM(IOP,I,J)
	  ENDDO
	ENDDO


	CALL INVMAT(BOX,INVBOX)
	CALL MULTMAT(EPS,INVBOX,R1)
	CALL MULTMAT(BOX,R1,R)


c     determine quaternions from rotational matrix
         
          aq = R(1,1)+R(2,2)
          bq = R(2,1)-R(1,2)
          cq = R(3,2)-R(2,3)
          dq = R(2,1)+R(1,2)
          eq = R(3,1)+R(1,3)
          fq = R(3,2)+R(2,3)
          gq = R(3,1)-R(1,3)
          hq = R(1,1)-R(2,2)
         
          q0 = 0.5d0*sqrt(aq+sqrt(aq*aq+bq*bq))
         
          if(q0.gt.1.d-4) then
          
            q1 = -0.25*cq/q0
            q2 =  0.25*gq/q0
            q3 = -0.25*bq/q0
          
          else
          
            q1 = 0.5*sqrt(hq+sqrt(hq*hq+dq*dq))
          
            if(q1.gt.1.d-4) then
          
              q2 = 0.25d0*dq/q1
              q3 = 0.25d0*eq/q1
          
            else
          
              q2 = 0.5d0*sqrt(-hq+sqrt(hq*hq+dq*dq))
          
              if(q2.gt.1.d-4) then
          
                q3 = 0.25*fq/q2
          
              else
          
                q3 = 1.
          
              endif
          
            endif
          
          endif
c
c     normalise quaternions
         
          rnorm = 1./sqrt(q0**2+q1**2+q2**2+q3**2)
          q0 =rnorm*q0
          q1 =rnorm*q1
          q2 =rnorm*q2
          q3 =rnorm*q3

      return
      end

C**************************************************

      SUBROUTINE INVMAT(A,B)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3),B(3,3),R(3,3)
      R(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3)
      R(2,1)=-A(2,1)*A(3,3)+A(3,1)*A(2,3)
      R(3,1)=A(2,1)*A(3,2)-A(3,1)*A(2,2)
      R(1,2)=-A(1,2)*A(3,3)+A(3,2)*A(1,3)
      R(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3)
      R(3,2)=-A(1,1)*A(3,2)+A(3,1)*A(1,2)
      R(1,3)=A(1,2)*A(2,3)-A(2,2)*A(1,3)
      R(2,3)=-A(1,1)*A(2,3)+A(2,1)*A(1,3)
      R(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2)
      D=A(1,1)*R(1,1)+A(1,2)*R(2,1)+A(1,3)*R(3,1)
      IF(ABS(D).LT.0.0001)
     *STOP'***Error INVMAT: wanishing matrix determinant'
      DO I=1,3
        DO J=1,3
          B(I,J)=R(I,J)/D
	ENDDO
      ENDDO
      RETURN
      END

C**************************************************
                                                                                                                                       
      SUBROUTINE MULTMATVEC(A,X_IN,Y_IN,Z_IN,X_OUT,Y_OUT,Z_OUT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3),B(3),C(3)
                                                                                                                                       
C       (X_OUT,Y_OUT,Z_OUT) = A.(X_IN,Y_IN,Z_IN)
                          
	X_OUT = A(1,1)*X_IN + A(1,2)*Y_IN + A(1,3)*Z_IN
	Y_OUT = A(2,1)*X_IN + A(2,2)*Y_IN + A(2,3)*Z_IN
	Z_OUT = A(3,1)*X_IN + A(3,2)*Y_IN + A(3,3)*Z_IN

      RETURN
      END
                                                                                                                                       
C**************************************************

C**************************************************

      SUBROUTINE MULTMAT(A,B,C)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3),B(3,3),C(3,3)

C 	C= A.B

	DO I=1,3
	  DO J=1,3
	    S = 0.0
	    DO K=1,3
	      S = S + A(I,K)*B(K,J) 
	    ENDDO
	    C(I,J) = S
	  ENDDO
	ENDDO

      RETURN
      END

C**************************************************

      SUBROUTINE MULTMAT_VEC(A,V,U)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(3,3),V(3),U(3)
C	U = A*V

        DO I=1,3
	  S = 0.0
          DO K=1,3
            S = S + A(I,K)*V(K) 
	  ENDDO
	  U(I) = S
        ENDDO

      RETURN
      END


C**************************************************

	SUBROUTINE MK_DLPOLY(MX,MY,MZ,NSPEC,NOP,NMOL,N_INTERACTION,
     $  OUT_DLPOLY_CONFIG,OUT_SYS_PDB,OUT_DLPOLY_FIELD,
     $  COMPOUND_NAME,MOLECULE_NAME,POTENTIEL_NAME,
     $  BOX,BOX_DIL,DILATATION,X,Y,Z,TYPE,NATMS,
     $  CHEMICAL,SITE_SPEC,
     $  MASS,CHARGE,INDEX_NB,POT_NB,TSYM,CHEM_TYPE,N_RIGID,
     $  N_CONST,N_BONDS,N_ANGLE,N_HARMD,N_TORSION,N_IMPRO,
     $                       LISTE_ATOM_RIGID,
     $                       INDEX_CONST,POT_CONST,
     $                       INDEX_ANGLE,POT_ANGLE,
     $                       INDEX_BONDS,POT_BONDS,
     $                       INDEX_TORSION,POT_TORSION,
     $                       INDEX_IMPRO,POT_IMPRO,
     $                       INDEX_HARMD,POT_HARMD,FF,
     $                       DLPOLY_TYPE,bigbox)

	IMPLICIT DOUBLE PRECISION(A-H,O-Z)

        PARAMETER(NATOM_MAX=100)
        PARAMETER(NSPEC_MAX=9)
        PARAMETER(NOP_MAX=108)
        PARAMETER(N=999999)
        PARAMETER(LEVCFG=0)

	CHARACTER*40 OUT_DLPOLY_CONFIG,OUT_DLPOLY_FIELD
	CHARACTER*40 OUT_SYS_PDB
        CHARACTER*60 COMPOUND_NAME
        CHARACTER*11 MOLECULE_NAME(NSPEC_MAX)
	CHARACTER*13 POTENTIEL_NAME
        DIMENSION BOX(3,3)
        DIMENSION BOX_DIL(3,3)
        DIMENSION X(NATOM_MAX,NSPEC_MAX)
        DIMENSION Y(NATOM_MAX,NSPEC_MAX)
        DIMENSION Z(NATOM_MAX,NSPEC_MAX)
        INTEGER TYPE(NATOM_MAX,NSPEC_MAX)
        DIMENSION NATMS(NATOM_MAX)
        CHARACTER*7 CHEMICAL(NATOM_MAX,NSPEC_MAX)
	INTEGER SITE_SPEC(NOP_MAX)
        REAL MASS(NATOM_MAX,NSPEC_MAX)
        DIMENSION CHARGE(NATOM_MAX,NSPEC_MAX)
        DIMENSION INDEX_NB(2,N),POT_NB(3,N)
        DIMENSION TSYM(NOP_MAX,3,4)
	CHARACTER*3 CHEM_TYPE(N)
        CHARACTER*7 DLPOLY_TYPE(N)
	DIMENSION LISTE_ATOM_RIGID(NATOM_MAX,NSPEC_MAX)
        DIMENSION INDEX_CONST(2,N,NSPEC_MAX),POT_CONST(1,N,NSPEC_MAX)
        DIMENSION INDEX_ANGLE(3,N,NSPEC_MAX),POT_ANGLE(2,N,NSPEC_MAX)
        DIMENSION INDEX_BONDS(2,N,NSPEC_MAX),POT_BONDS(2,N,NSPEC_MAX)
        DIMENSION INDEX_TORSION(4,N,NSPEC_MAX)
        DIMENSION POT_TORSION(5,N,NSPEC_MAX)
        DIMENSION INDEX_HARMD(4,N,NSPEC_MAX),POT_HARMD(2,N,NSPEC_MAX)
        DIMENSION INDEX_IMPRO(4,N,NSPEC_MAX),POT_IMPRO(2,N,NSPEC_MAX)

	CHARACTER*5 FF

	DIMENSION BBOX(3,3)
	DIMENSION RR(3,3)
	INTEGER   ZSPEC(NSPEC_MAX)
        REAL MK
        INTEGER TK
	LOGICAL OK
	LOGICAL OK2,bigbox

        DIMENSION  N_RIGID(NSPEC_MAX),N_CONST(NSPEC_MAX)
        DIMENSION  N_ANGLE(NSPEC_MAX),N_TORSION(NSPEC_MAX)
        DIMENSION  N_HARMD(NSPEC_MAX),N_BONDS(NSPEC_MAX)
        DIMENSION  N_IMPRO(NSPEC_MAX)



        DO I=1,NSPEC_MAX
          ZSPEC(I) = 0
        ENDDO
        DO I=1,NOP
          ZSPEC(SITE_SPEC(I)) = ZSPEC(SITE_SPEC(I)) + 1
        ENDDO


        BX = REAL(MX)
        BY = REAL(MY)
        BZ = REAL(MZ)

        DO I=1,3
            BBOX(I,1) = BOX(1,I)
            BBOX(I,2) = BOX(2,I)
            BBOX(I,3) = BOX(3,I)
        ENDDO
        DO I=1,3
            BBOX(1,I) = BBOX(1,I)*BX
            BBOX(2,I) = BBOX(2,I)*BY
            BBOX(3,I) = BBOX(3,I)*BZ
        ENDDO

        IF( (BOX(1,2) .EQ. 0.0) .AND.
     $      (BOX(1,3) .EQ. 0.0) .AND.
     $      (BOX(2,3) .EQ. 0.0)  ) THEN
          IF(BOX(1,1) .EQ. BOX(2,2) 
     $       .AND. BOX(1,1) .EQ. BOX(3,3)) THEN
            IMCON=1
          ELSE
            IMCON=2
          ENDIF
        ELSE
          IMCON=3
        ENDIF

C       MK DL_POLY CONFIG

        OPEN(UNIT=2,FILE=OUT_DLPOLY_CONFIG,STATUS='UNKNOWN')
        WRITE(2,'(a40)') COMPOUND_NAME
        WRITE(2,'(2I10)')  LEVCFG,IMCON

	if(bigbox) then
          WRITE(2,'(3F20.10)') 50.0,0.0,0.0
          WRITE(2,'(3F20.10)') 0.0,50.0,0.0
          WRITE(2,'(3F20.10)') 0.0,0.0,50.0
	else
          WRITE(2,'(3F20.10)') BBOX(1,1),BBOX(1,2),BBOX(1,3)
          WRITE(2,'(3F20.10)') BBOX(2,1),BBOX(2,2),BBOX(2,3)
          WRITE(2,'(3F20.10)') BBOX(3,1),BBOX(3,2),BBOX(3,3)
	endif


        INDEX=0
	N_MOLECULE=0
        OK = .TRUE.
        DO ISPEC=1,NSPEC
	  DO IOP=1,NOP
	    IF(SITE_SPEC(IOP) .EQ. ISPEC) THEN 
	      DO IX=1,MX
	        DO IY=1,MY
	          DO IZ=1,MZ

                    XMAILLE = REAL(IX-1)
                    YMAILLE = REAL(IY-1)
                    ZMAILLE = REAL(IZ-1)

	             N_MOLECULE=N_MOLECULE+1

	             IF(OK) THEN 
                     DO K=1,NATMS(ISPEC)
                       TK = TYPE(K,ISPEC)
                       XK = X(K,ISPEC)
                       YK = Y(K,ISPEC)
                       ZK = Z(K,ISPEC)

                         XTAL = XK*TSYM(IOP,1,1) + 
     $                          YK*TSYM(IOP,1,2) +
     $                          ZK*TSYM(IOP,1,3) + 
     $                          TSYM(IOP,1,4) + XMAILLE
                         YTAL = XK*TSYM(IOP,2,1) + 
     $                          YK*TSYM(IOP,2,2) +
     $                          ZK*TSYM(IOP,2,3) + 
     $                          TSYM(IOP,2,4) + YMAILLE
                         ZTAL = XK*TSYM(IOP,3,1) + 
     $                          YK*TSYM(IOP,3,2) +
     $                          ZK*TSYM(IOP,3,3) + 
     $                          TSYM(IOP,3,4) + ZMAILLE


		       CALL MULTMATVEC(BOX,XTAL,YTAL,ZTAL,XO,YO,ZO)

                       INDEX = INDEX + 1


                       WRITE(2,'(a,I10)') CHEMICAL(K,ISPEC),INDEX
                       WRITE(2,'(3F20.10)') XO,YO,ZO


                     ENDDO ! NATMS
	             ENDIF ! OK

	          ENDDO ! MZ 
	        ENDDO ! MY
	      ENDDO ! MX
	    ENDIF ! SITE_SPEC
          ENDDO ! NOP
        ENDDO ! NSPEC

        CLOSE(2)
        PRINT *,'mk Config  ok'


C       MK DL_POLY FIELD

        OPEN(UNIT=2,FILE=OUT_DLPOLY_FIELD,STATUS='UNKNOWN')
        WRITE(2,'(a40)') COMPOUND_NAME
        WRITE(2,'(a)') 'UNITS kcal'
        WRITE(2,'(a9,1X,I5)') 'MOLECULES',NSPEC
        DO ISPEC=1,NSPEC
          WRITE(2,'(a11)') MOLECULE_NAME(ISPEC)
C       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C	ATTENTION NE MARCHE QU'AVEC 1 TYPE DE MOLECULES !!!
C       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          WRITE(2,'(a7,I5)') 'nummols',NMOL*ZSPEC(ISPEC)
          WRITE(2,'(a5,I5)') 'atoms',NATMS(ISPEC)
          DO K=1,NATMS(ISPEC)
            TK = TYPE(K,ISPEC)
            WRITE(2,'(A8,2F12.6,2I5)')
     *     CHEMICAL(K,ISPEC),MASS(K,ISPEC),CHARGE(K,ISPEC),1,0
          ENDDO
          IF(N_RIGID(ISPEC) .GT. 0) THEN
            WRITE(2,'(a5,I5)') 'rigid',1
            WRITE(2,'(I5,$)') N_RIGID(ISPEC)
            K = 1
            DO I=1,N_RIGID(ISPEC)
              WRITE(2,'(I5,$)') LISTE_ATOM_RIGID(I,ISPEC)
              K = K + 1
              IF(MOD(K,16) .EQ. 0) WRITE(2,*)
            ENDDO
            WRITE(2,*)
          ENDIF

          IF(N_CONST(ISPEC) .GT. 0) THEN
            WRITE(2,'(a11,1x,I5)') 'constraints',N_CONST(ISPEC)
            DO I=1,N_CONST(ISPEC)
               WRITE(2,'(2I5,F12.6)') 
     $        (INDEX_CONST(K,I,ISPEC),K=1,2),
     $         POT_CONST(1,I,ISPEC)
            ENDDO
          ENDIF


          IF(N_BONDS(ISPEC) .GT. 0) THEN
            WRITE(2,'(a6,I5)') 'bonds',N_BONDS(ISPEC)
            DO I=1,N_BONDS(ISPEC)
               WRITE(2,'(a4,2I5,2F12.6)') 'harm-hrm',
     $        (INDEX_BONDS(K,I,ISPEC),K=1,2),
     $         POT_BONDS(1,I,ISPEC),POT_BONDS(2,I,ISPEC)
            ENDDO
          ENDIF

          IF(N_ANGLE(ISPEC) .GT. 0) THEN
            WRITE(2,'(a6,I5)') 'angles',N_ANGLE(ISPEC)
            DO I=1,N_ANGLE(ISPEC)
               WRITE(2,'(a4,3I5,2F12.6)') 'harm',
     $        (INDEX_ANGLE(K,I,ISPEC),K=1,3),
     $         POT_ANGLE(1,I,ISPEC),POT_ANGLE(2,I,ISPEC)
            ENDDO
          ENDIF

          IF(N_TORSION(ISPEC) .GT. 0 .OR. N_IMPRO(ISPEC) .GT. 0 ) THEN
            WRITE(2,'(a9,I5)') 'dihedrals',
     $      (N_TORSION(ISPEC)+N_IMPRO(ISPEC))

            IF(N_TORSION(ISPEC) .GT. 0) THEN
            DO I=1,N_TORSION(ISPEC)
               WRITE(2,'(a3,4I5,5F12.6)') 'cos',
     $        (INDEX_TORSION(K,I,ISPEC),K=1,4),
     $        POT_TORSION(1,I,ISPEC),POT_TORSION(3,I,ISPEC)
     $        ,POT_TORSION(2,I,ISPEC),
     $       POT_TORSION(4,I,ISPEC),POT_TORSION(5,I,ISPEC)
            ENDDO
            ENDIF

            IF(N_IMPRO(ISPEC) .GT. 0) THEN
            DO I=1,N_IMPRO(ISPEC)
               WRITE(2,'(a3,4I5,5F12.6)') 'cos',
     $        (INDEX_IMPRO(K,I,ISPEC),K=1,4),
     $        0.550000,180.000000,2.000000,0.0,0.0
            ENDDO

            ENDIF


          ENDIF
          WRITE(2,'(a6)') 'finish'
        ENDDO
        WRITE(2,'(a3,I5)') 'VDW',N_INTERACTION

	IF (FF .EQ. 'GAFF_') THEN
          DO I=1,N_INTERACTION
            K = INDEX_NB(1,I)
            L = INDEX_NB(2,I)
            A = POT_NB(1,I)
            B = POT_NB(2,I)
            WRITE(2,'(a7,1x,a7,1x,a4,X,F12.3,F12.6)')
     $      CHEM_TYPE(K),CHEM_TYPE(L),'lj',B,A
          ENDDO
        ELSE IF(FF .EQ. 'OPLS_') THEN
          DO I=1,N_INTERACTION
            K = INDEX_NB(1,I)
            L = INDEX_NB(2,I)
            A = POT_NB(1,I)
            B = POT_NB(2,I)
            WRITE(2,'(a7,1x,a7,1x,a4,X,F12.6,F12.6)')
     $      DLPOLY_TYPE(K),DLPOLY_TYPE(L),'lj',B,A
          ENDDO
	ELSE
	  PRINT *,'Unkwown FF in MK_DLPOLY'
	  STOP
	ENDIF ! FF

        WRITE(2,'(a5)') 'CLOSE'
        CLOSE(2)
        PRINT *,'mk Field  ok'

	RETURN
	END

C**************************************************

C       BONDS
        FUNCTION BONDS(X1,Y1,Z1,X2,Y2,Z2,BOX)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION BOX(3,3)
        PI = ACOS(-1.)

        DXA = X1 - X2
        DYA = Y1 - Y2
        DZA = Z1 - Z2

        DXA = DXA -ANINT(DXA)
        DYA = DYA -ANINT(DYA)
        DZA = DZA -ANINT(DZA)

	CALL MULTMATVEC(BOX,DXA,DYA,DZA,DX1,DY1,DZ1)

        BONDS  = SQRT(DX1**2+DY1**2+DZ1**2)

        END

C*************************************************************

C       ANGLE
        FUNCTION ANGLE(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,BOX)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION BOX(3,3)
        PI = ACOS(-1.)
        RD = 180./PI

        DXA = X1 - X2
        DYA = Y1 - Y2
        DZA = Z1 - Z2

        DXB = X3 - X2
        DYB = Y3 - Y2
        DZB = Z3 - Z2

        DXA = DXA -ANINT(DXA)
        DYA = DYA -ANINT(DYA)
        DZA = DZA -ANINT(DZA)

        DXB = DXB -ANINT(DXB)
        DYB = DYB -ANINT(DYB)
        DZB = DZB -ANINT(DZB)

	CALL MULTMATVEC(BOX,DXA,DYA,DZA,DX1,DY1,DZ1)
	CALL MULTMATVEC(BOX,DXB,DYB,DZB,DX2,DY2,DZ2)

        R1  = SQRT(DX1**2+DY1**2+DZ1**2)
        DX1 = DX1/R1
        DY1 = DY1/R1
        DZ1 = DZ1/R1
        R2  = SQRT(DX2**2+DY2**2+DZ2**2)
        DX2 = DX2/R2
        DY2 = DY2/R2
        DZ2 = DZ2/R2

        ANGLE = RD*ACOS(DX1*DX2 + DY1*DY2 + DZ1*DZ2)

        END

C**************************************************

C       DIHED
        FUNCTION DIHED(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,BOX)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION BOX(3,3)
        PI = ACOS(-1.)
        RD = 180./PI

        DXA = X2 - X1
        DYA = Y2 - Y1
        DZA = Z2 - Z1

        DXB = X3 - X2
        DYB = Y3 - Y2
        DZB = Z3 - Z2

        DXC = X4 - X3
        DYC = Y4 - Y3
        DZC = Z4 - Z3

        DXA = DXA -ANINT(DXA)
        DYA = DYA -ANINT(DYA)
        DZA = DZA -ANINT(DZA)

        DXB = DXB -ANINT(DXB)
        DYB = DYB -ANINT(DYB)
        DZB = DZB -ANINT(DZB)

        DXC = DXC -ANINT(DXC)
        DYC = DYC -ANINT(DYC)
        DZC = DZC -ANINT(DZC)

	CALL MULTMATVEC(BOX,DXA,DYA,DZA,XA,YA,ZA)
	CALL MULTMATVEC(BOX,DXB,DYB,DZB,XB,YB,ZB)
	CALL MULTMATVEC(BOX,DXC,DYC,DZC,XC,YC,ZC)

        DX1 = YA*ZB - YB*ZA
        DY1 = ZA*XB - ZB*XA
        DZ1 = XA*YB - YA*XB

        DX2 = YB*ZC - YC*ZB
        DY2 = ZB*XC - ZC*XB
        DZ2 = XB*YC - YB*XC

        R1  = SQRT(DX1**2+DY1**2+DZ1**2)
        DX1 = DX1/R1
        DY1 = DY1/R1
        DZ1 = DZ1/R1

        R2  = SQRT(DX2**2+DY2**2+DZ2**2)
        DX2 = DX2/R2
        DY2 = DY2/R2
        DZ2 = DZ2/R2 

        DIHED = RD*ACOS(DX1*DX2 + DY1*DY2 + DZ1*DZ2)

        END


	SUBROUTINE ORTH2CRYS(X,Y,Z,BOX)
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	DIMENSION BOX(3,3)

	CALL MULTMATVEC(BOX,X,Y,Z,XN,YN,ZN)
	X = XN
	Y = YN
	Z = ZN

	RETURN
	END

	SUBROUTINE CAL_DIPOLE(NATOM_MAX,NSPEC_MAX,NSPEC,
     $  NATMS,X,Y,Z,BOX,CHARGE,MOLECULE_NAME)
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	PARAMETER(UNIT_DEBYE = 4.802)
        DIMENSION X(NATOM_MAX,NSPEC_MAX)
        DIMENSION Y(NATOM_MAX,NSPEC_MAX)
        DIMENSION Z(NATOM_MAX,NSPEC_MAX)
	DIMENSION CHARGE(NATOM_MAX,NSPEC_MAX)
	DIMENSION NATMS(NATOM_MAX)
	DIMENSION BOX(3,3)
	CHARACTER*11 MOLECULE_NAME(NSPEC_MAX)

        DO ISPEC=1,NSPEC
	  DX = 0.
	  DY = 0.
	  DZ = 0.
	  PRINT *,'spec = ',MOLECULE_NAME(ISPEC)
          DO K=1,NATMS(ISPEC)
            XTAL = X(K,ISPEC)
            YTAL = Y(K,ISPEC)
            ZTAL = Z(K,ISPEC)
	    CALL MULTMATVEC(BOX,XTAL,YTAL,ZTAL,XO,YO,ZO)
	    DX = DX + XO*CHARGE(K,ISPEC)
	    DY = DY + YO*CHARGE(K,ISPEC)
	    DZ = DZ + ZO*CHARGE(K,ISPEC)
	  ENDDO
	  DX = DX*UNIT_DEBYE
	  DY = DY*UNIT_DEBYE
	  DZ = DZ*UNIT_DEBYE
	  D2 = SQRT(DX**2 + DY**2 + DZ**2)
	  PRINT '(a,F10.5)','mu   = ',D2
	  PRINT '(a,F10.5)','mu_x = ',DX   
	  PRINT '(a,F10.5)','mu_y = ',DY 
	  PRINT '(a,F10.5)','mu_z = ',D2
	  PRINT *,' '
	ENDDO

	RETURN
	END


	SUBROUTINE CHECK_IF_BOND_EXIST(I1,I2,ISPEC,N,NSPEC_MAX,
     $  INDEX_BONDS,INDEX_CONST,N_BONDS,N_CONST)
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 	DIMENSION INDEX_BONDS(2,N,NSPEC_MAX)
 	DIMENSION INDEX_CONST(2,N,NSPEC_MAX)
	DIMENSION N_BONDS(NSPEC_MAX)
	DIMENSION N_CONST(NSPEC_MAX)
	LOGICAL BOND_EXIST



        BOND_EXIST = .FALSE.
        IF(N_BONDS(ISPEC) .GT. 0) THEN
           DO JJ=1,N_BONDS(ISPEC)
             II1 = INDEX_BONDS(1,JJ,ISPEC)
             II2 = INDEX_BONDS(2,JJ,ISPEC)
             IF( ((I1 .EQ. II1) .AND. (I2 .EQ. II2)) .OR.
     $           ((I1 .EQ. II2) .AND. (I2 .EQ. II1)) ) THEN
                BOND_EXIST = .TRUE.
             ENDIF
           ENDDO
         ELSEIF(N_CONST(ISPEC) .GT. 0) THEN
           DO JJ=1,N_CONST(ISPEC)
             II1 = INDEX_CONST(1,JJ,ISPEC)
             II2 = INDEX_CONST(2,JJ,ISPEC)
             IF( ((I1 .EQ. II1) .AND. (I2 .EQ. II2)) .OR.
     $           ((I1 .EQ. II2) .AND. (I2 .EQ. II1)) ) THEN
               BOND_EXIST = .TRUE.
             ENDIF
           ENDDO
	 ELSE
	   PRINT *,'no bond defined !!!'
	   STOP
         ENDIF

	IF(.NOT. BOND_EXIST) THEN
	  PRINT '(a,i5,a,i5)','bond : ',I1,' - ',I2
	  PRINT *,'bond does not exist !!!'
	  STOP
	ENDIF
	
	RETURN
	END


        SUBROUTINE GENERATE_ANGLE_LISTE(ISPEC,N,NSPEC_MAX,
     $  INDEX_B,NB,NANGLE,INDEX_A,NATOM)
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION INDEX_B(2,N,NSPEC_MAX)
	DIMENSION INDEX_A(3,N,NSPEC_MAX)
	DIMENSION LIST(NATOM)

C	FILL ANGLE LIST  (X - IATOM - Y)
	J = 0
	DO IATOM = 1,NATOM
	  K = 0
C FOUND ALL ATOMS BOUND TO IATOM 
	  DO IB=1,NB
	    I1 = INDEX_B(1,IB,ISPEC)
	    I2 = INDEX_B(2,IB,ISPEC)
	    IF(IATOM .EQ. I1) THEN 
	      K = K + 1
	      LIST(K) = I2
	    ELSEIF (IATOM .EQ. I2 ) THEN
	      K = K + 1
	      LIST(K) = I1
	    ENDIF
	  ENDDO  
	  IF(K .GE. 2) THEN 
	    DO IK1=1,K-1
	      DO IK2=IK1+1,K
	        J = J + 1
	        INDEX_A(1,J,ISPEC) = LIST(IK1)
	        INDEX_A(2,J,ISPEC) = IATOM
	        INDEX_A(3,J,ISPEC) = LIST(IK2)
	      ENDDO
	    ENDDO
	  ENDIF
	ENDDO 

	NANGLE = J


	RETURN
	END


        SUBROUTINE GENERATE_TORSI_LISTE(ISPEC,N,NSPEC_MAX,
     $  INDEX_B,NB,NTORSI,INDEX_T,NATOM)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION INDEX_B(2,N,NSPEC_MAX)
        DIMENSION INDEX_T(4,N,NSPEC_MAX)
        DIMENSION LIST1(NATOM)
        DIMENSION LIST2(NATOM)

C       FILL TORSION LIST  (X - BONDS - Y)

        J = 0
	DO IB=1,NB
          IB1 = INDEX_B(1,IB,ISPEC)
          IB2 = INDEX_B(2,IB,ISPEC)
C FOUND ALL ATOMS BOUND TO IB1 EXECEPT IB2
	  K = 0
          DO I=1,NB
            I1 = INDEX_B(1,I,ISPEC)
            I2 = INDEX_B(2,I,ISPEC)
            IF(IB1 .EQ. I1 .AND. I2 .NE. IB2) THEN
              K= K+ 1
              LIST1(K) = I2
            ELSEIF (IB1 .EQ. I2 .AND. I1 .NE. IB2) THEN
              K = K + 1
              LIST1(K) = I1
            ENDIF
          ENDDO
	  K1 = K
C FOUND ALL ATOMS BOUND TO IB2 EXECEPT IB1
          K = 0
          DO I=1,NB
            I1 = INDEX_B(1,I,ISPEC)
            I2 = INDEX_B(2,I,ISPEC)
            IF(IB2 .EQ. I1 .AND. I2 .NE. IB1) THEN
              K= K+ 1
              LIST2(K) = I2
            ELSEIF (IB2 .EQ. I2 .AND. I1 .NE. IB1) THEN
              K = K + 1
              LIST2(K) = I1
            ENDIF
          ENDDO
          K2 = K

          IF( (K1 .GE. 1) .AND. (K2 .GE. 1) ) THEN
            DO IK1=1,K1
              DO IK2=1,K2
                J = J + 1
                INDEX_T(1,J,ISPEC) = LIST1(IK1)
                INDEX_T(2,J,ISPEC) = IB1
                INDEX_T(3,J,ISPEC) = IB2
                INDEX_T(4,J,ISPEC) = LIST2(IK2)
              ENDDO
            ENDDO
          ENDIF

	ENDDO
        NTORSI = J



        RETURN
        END



C**************************************************
C	FILL TAB
C**************************************************

	SUBROUTINE FF_PARAM(FF,INTRA_PROP,C1,C2,C3,C4,
     $      IP,III,PARAM,longchain,TYPE_TORSION)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
	PARAMETER(N=999999)
	DIMENSION X(N)
	DIMENSION PARAM(N)
	CHARACTER*1 DUM
	CHARACTER*3 C1,C2,C3,C4
	CHARACTER*3 CR1,CR2,CR3,CR4
	CHARACTER*3 SEC1(n),SEC2(n),SEC3(n),SEC4(n)
	CHARACTER*5 INTRA_PROP,FF
	CHARACTER*35 FF_FILE
	CHARACTER*2 ATOM
	CHARACTER*80 HEADER
	CHARACTER*20  shortchain
	CHARACTER*20  chainmatch(n)
	CHARACTER*20  TYPE_TORSION(n)
	CHARACTER*130 longchain
	LOGICAL OK,OK2(n)
        REAL*8 MASS_ATOM
        integer listmatch(N),iposchain(4)
	OK = .TRUE.

C       FORCE_FIELD
        FF_file(1:20) = '/mnt/c/bin/FF'
        FF_file(21:25) = FF
        FF_file(26:26) = '/'
        FF_file(27:31) = INTRA_PROP
        FF_file(32:35) = '.tab'

        do i=1,n
          PARAM(i) = 0.0
        enddo 

C*********************
C************* BONDS
C*********************
	IF(INTRA_PROP .EQ. 'bonds') THEN
	  OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
	  READ(11,'(a)') HEADER
          DO I=1,N
            READ(11,*,end=100) CR1,CR2,(X(J),J=1,2)
c           READ(11,'(a3,a,a3,2F10.5)',end=100) CR1,DUM,CR2,(X(J),J=1,2)
c               PRINT *,"I = ",I
c               PRINT *,"CR1",CR1
c   PRINT *,"DUM",DUM
c       PRINT *,"CR2",CR2
c       PRINT *,"X1 ",X(1)
c       PRINT *,"X2 ",X(2)
            IF((CR1 .EQ. C1 .AND. CR2 .EQ. C2) .OR.
     $         (CR1 .EQ. C2 .AND. CR2 .EQ. C1) 
     $         .AND. OK ) THEN 
              PARAM(1) = X(1)
              PARAM(2) = X(2)
	      OK = .FALSE.
            ENDIF
	  ENDDO
100	  CLOSE(11)
          IF(OK) THEN 
            PRINT *,'Unknown bonds interaction = ',
     $       C1,'-',C2
	    PRINT *,'Upgrade bonds.tab !!!'
            STOP
          ENDIF


C*********************
C************* ANGLE
C*********************
	ELSEIF(INTRA_PROP .EQ. 'angle') THEN
	  OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
	  READ(11,'(a)') HEADER
          DO I=1,N
            READ(11,*,end=200) CR1,CR2,CR3,(X(J),J=1,2)
c           READ(11,'(a3,a,a3,a,a3,2f10.5)',end=200) 
c    $        CR1,DUM,CR2,DUM,CR3,(X(J),J=1,2)
            IF((CR1 .EQ. C1 .AND. CR2 .EQ. C2 
     $                      .AND. CR3 .EQ. C3) .OR.
     $         (CR1 .EQ. C3 .AND. CR2 .EQ. C2 
     $                      .AND. CR3 .EQ. C1)
     $         .AND. OK ) THEN
              PARAM(1) = X(1)
              PARAM(2) = X(2)
              OK = .FALSE.
            ENDIF
          ENDDO
200	  CLOSE(11)
          IF(OK) THEN 
            PRINT *,'Unknown angle interaction = ',
     $      C1,'-',C2,'-',C3
	    PRINT *,'Upgrade angle.tab !!!'
            STOP
          ENDIF


C*********************
C************* TORSION 
C*********************

	ELSEIF(INTRA_PROP .EQ. 'torsi') THEN


        n_type_torsion = III
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       search for exact match
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


	    match_ok = 0
            kk = 0
	    OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
	    READ(11,'(a)') HEADER
            DO I=1,N
              READ(11,*,end=300) CR1,CR2,CR3,CR4,X(1),X(2),X(3),
     $         X(4),shortchain(1:19)
              shortchain(20:20) = "#"


	      IF( ( (CR1 .EQ. C1) .AND. 
     $              (CR2 .EQ. C2) .AND.
     $              (CR3 .EQ. C3) .AND. 
     $              (CR4 .EQ. C4) ) .OR.
     $            ( (CR4 .EQ. C1) .AND.
     $              (CR3 .EQ. C2) .AND.
     $              (CR2 .EQ. C3) .AND.
     $              (CR1 .EQ. C4) ) ) THEN 


c               print *,C1,C2,C3,C4
c               print *,CR1,CR2,CR3,CR4,shortchain
c               print *,' '
	        match_ok = match_ok + 1
                PARAM(kk+1) = X(1)
                PARAM(kk+2) = X(2)
                PARAM(kk+3) = X(3)
                kk = kk + 3 
                listmatch(match_ok) = I
                chainmatch(match_ok) = shortchain
                OK = .FALSE.

	      ENDIF
            ENDDO
300	    CLOSE(11)

            IF(.NOT. OK) THEN 

C       eliminate multiple solution using torsion_type list

	    if(match_ok .gt. 1) then
              match_ok2 = 0
              DO kk = 1,match_ok
                OK2(kk) = .false. 

                DO j=1,n_type_torsion
                 if(chainmatch(kk) .eq. type_torsion(j) ) 
     $            ok2(kk) = .true.
                ENDDO

                if(ok2(kk)) then 
                  match_ok2 = match_ok2 + 1
                  kk2 = kk
                endif 
              ENDDO
            else
                match_ok2 =  match_ok
                kk2 = 1
            endif 

	    if(match_ok2 .gt. 1) then
	      print *,'exact match solution '
              print *,'!!!!!!! several solutions for : ',
     $        C1,"-",C2,"-",C3,"-",C4,' : ',match_ok2
              do kk=1,match_ok
                 if(ok2(kk)) then 
                   P1 = PARAM(1 + 3*(kk-1))
                   P2 = PARAM(2 + 3*(kk-1))
                   P3 = PARAM(3 + 3*(kk-1))
                   print '(2i5,3F8.3,2a)',kk,listmatch(kk),
     $             P1,P2,P3,' : ',chainmatch(kk)
                 endif 
              enddo 
            elseif(match_ok2 .eq. 1) then 
	      print *,'exact match solution '
	      print *,'one solution for : ',
     $        C1,"-",C2,"-",C3,"-",C4,'  '

                   P1 = PARAM(1 + 3*(kk2-1))
                   P2 = PARAM(2 + 3*(kk2-1))
                   P3 = PARAM(3 + 3*(kk2-1))
                print '(i5,3F8.3,2a)',listmatch(kk2),
     $             P1,P2,P3,' : ',chainmatch(kk2)
                PARAM(1) = P1
                PARAM(2) = P2
                PARAM(3) = P3
                chainmatch(1) = chainmatch(kk2)
                listmatch(1) = listmatch(kk2)
            else
              print *,'!! add torsion type for'
	      print *,
     $        C1,"-",C2,"-",C3,"-",C4,'  '
              do kk=1,match_ok
                 print '(a20)',chainmatch(kk)
              enddo
	    endif 
              print *,' '
           ENDIF 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C       search for partial solution C* O* * 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

          IF(OK) THEN 
C       search for other match C* O* ... 
               print *,'search for partial solutions'
               print '(4a)',C1(1:3),C2(1:3),C3(1:3),C4(1:3)
        
            match_ok = 0
            kk = 0
            OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
            READ(11,'(a)') HEADER
            DO I=1,N
              READ(11,*,end=310) CR1,CR2,CR3,CR4,X(1),X(2),X(3),
     $         X(4),shortchain(1:19)
              shortchain(20:20) = "#"
               IF(
     $            (
     $             (  
     $               (CR1(1:1) .EQ. '*') .OR. 
     $               (CR1 .EQ. C1) .OR. 
     $               ((CR1(1:1) .EQ. C1(1:1)) 
     $               .AND. (CR1(2:2) .EQ. '*' ))  
     $             ) 
     $             .AND. 
     $             (  
     $               (CR2(1:1) .EQ. '*') .OR. 
     $               (CR2 .EQ. C2) .OR. 
     $               ((CR2(1:1) .EQ. C2(1:1)) 
     $               .AND. (CR2(2:2) .EQ. '*' ))  
     $             ) 
     $             .AND.
     $             (  
     $               (CR3(1:1) .EQ. '*') .OR. 
     $               (CR3 .EQ. C3) .OR. 
     $               ((CR3(1:1) .EQ. C3(1:1)) 
     $               .AND. (CR3(2:2) .EQ. '*' ))  
     $             ) 
     $             .AND. 
     $             (  
     $               (CR4(1:1) .EQ. '*') .OR. 
     $               (CR4 .EQ. C4) .OR. 
     $               ((CR4(1:1) .EQ. C4(1:1)) 
     $               .AND. (CR4(2:2) .EQ. '*' ))  
     $             )
     $           ) .OR. 

     $            (
     $             (
     $               (CR4(1:1) .EQ. '*') .OR. 
     $               (CR4 .EQ. C1) .OR.
     $               ((CR4(1:1) .EQ. C1(1:1))
     $               .AND. (CR4(2:2) .EQ. '*' ))
     $             )
     $             .AND.
     $             (
     $               (CR3(1:1) .EQ. '*') .OR. 
     $               (CR3 .EQ. C2) .OR.
     $               ((CR3(1:1) .EQ. C2(1:1))
     $               .AND. (CR3(2:2) .EQ. '*' ))
     $             )
     $             .AND.
     $             (
     $               (CR2(1:1) .EQ. '*') .OR. 
     $               (CR2 .EQ. C3) .OR.
     $               ((CR2(1:1) .EQ. C3(1:1))
     $               .AND. (CR2(2:2) .EQ. '*' ))
     $             )
     $             .AND.
     $             (
     $               (CR1(1:1) .EQ. '*') .OR. 
     $               (CR1 .EQ. C4) .OR.
     $               ((CR1(1:1) .EQ. C4(1:1))
     $               .AND. (CR1(2:2) .EQ. '*' ))
     $             )
     $           ) 
     $            ) THEN 

                PARAM(kk+1) = X(1)
                PARAM(kk+2) = X(2)
                PARAM(kk+3) = X(3)
                kk = kk + 3
                match_ok = match_ok + 1
                listmatch(match_ok) = I
                chainmatch(match_ok) = shortchain
                SEC1(match_ok) = CR1
                SEC2(match_ok) = CR2
                SEC3(match_ok) = CR3
                SEC4(match_ok) = CR4
                OK = .FALSE.
              ENDIF
            ENDDO
310         CLOSE(11)


C       eliminate multiple solution using torsion_type list

            IF(.NOT. OK) THEN

            if(match_ok .gt. 1) then
              match_ok2 = 0
              DO kk = 1,match_ok
                OK2(kk) = .false.

                DO j=1,n_type_torsion
                 if(chainmatch(kk) .eq. type_torsion(j) )
     $            ok2(kk) = .true.
                   
                ENDDO

                if(ok2(kk)) then
                  match_ok2 = match_ok2 + 1
                  kk2 = kk
                endif
              ENDDO
            else
                match_ok2 =  match_ok
                kk2 = 1
            endif

            if(match_ok2 .gt. 1) then
              print *,'generic match *'
              print *,'!!!!!!! several solutions for : ',
     $        C1,"-",C2,"-",C3,"-",C4,' : ',match_ok2
              do kk=1,match_ok
                 if(ok2(kk)) then
                   P1 = PARAM(1 + 3*(kk-1))
                   P2 = PARAM(2 + 3*(kk-1))
                   P3 = PARAM(3 + 3*(kk-1))
                   print '(2i5,3F8.3,2a)',kk,listmatch(kk),
     $             P1,P2,P3,' : ',chainmatch(kk)
                 endif
              enddo
            elseif(match_ok2 .eq. 1) then
              print *,'generic match *'
              print *,'one solution for : ',
     $        C1,"-",C2,"-",C3,"-",C4,'  '

                   P1 = PARAM(1 + 3*(kk2-1))
                   P2 = PARAM(2 + 3*(kk2-1))
                   P3 = PARAM(3 + 3*(kk2-1))
                print '(i5,3F8.3,2a)',listmatch(kk2),
     $             P1,P2,P3,' : ',chainmatch(kk2)
                PARAM(1) = P1
                PARAM(2) = P2
                PARAM(3) = P3
                chainmatch(1) = chainmatch(kk2)
                listmatch(1) = listmatch(kk2)
            else
              print *,'!!!!!!!  add torsion type for'
              print *,
     $        C1,"-",C2,"-",C3,"-",C4,'  '
              do kk=1,match_ok
                 print '(a20)',chainmatch(kk)
              enddo
            endif
              print *,' '
           ENDIF

          ENDIF 

          IF(OK) THEN
            PRINT *,'!!!!!!  Unknown torsion interaction = ',
     $      C1,'-',C2,'-',C3,'-',C4 
            PRINT *,'Interaction potential set to zero  !!!'
            PRINT *,'Upgrade torsi.tab !!!'
            print *,' '
          ENDIF

          IP = match_ok2
C*********************
C************* MASSE & TYPE
C*********************

	ELSEIF(INTRA_PROP .EQ. 'masse') THEN

          IVDW = III
	  OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
	  READ(11,'(a)') HEADER
          DO I=1,N
            READ(11,*,end=500) itabvdw,CR1,ATOM
	    IF( (itabvdw .eq. ivdw) .AND. OK) THEN 
                PARAM(1) =  mass_atom(ATOM)
                NOK = I 
                OK = .false. 
	    ENDIF
          ENDDO
500	  CLOSE(11)

          IF(.NOT. OK) THEN 
            OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
            READ(11,'(a)') header
            DO I=1,NOK
              READ(11,'(a130)') longchain
            ENDDO
            CLOSE(11)
          ENDIF 

          IF(OK) THEN
            PRINT *,'Unknown chemical mass = ',C1
            PRINT *,'Upgrade masse.tab !!!'
            STOP
          ENDIF

C*********************
C************* VDW 
C*********************

        ELSEIF(INTRA_PROP .EQ. 'vdw1a') THEN
          IVDW = III
          OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
          READ(11,'(a)') HEADER
          DO I=1,N
            READ(11,*,end=600) itabvdw,CR1,X(1),X(2),X(3)
            IF( (itabvdw .eq. ivdw) .AND. OK) THEN 
                PARAM(1) = X(1)
                PARAM(2) = X(2)
                PARAM(3) = X(3)
                OK = .FALSE.
            ENDIF
          ENDDO
600       CLOSE(11)


C*********************
C************* Q charge 
C*********************

        ELSEIF(INTRA_PROP .EQ. 'vdw1a') THEN
          IVDW = III
          OPEN(UNIT=11,FILE=FF_file,STATUS='OLD')
          READ(11,'(a)') HEADER
          DO I=1,N
            READ(11,*,end=610) itabvdw,CR1,X(1)
            IF( (itabvdw .eq. ivdw) .AND. OK) THEN
                PARAM(1) = X(1)
                OK = .FALSE.
            ENDIF
          ENDDO
610       CLOSE(11)

          IF(OK) THEN
            PRINT *,'Unknown chemical FF = ',C1
            PRINT *,'Upgrade vdw1a.tab !!!'
            STOP
          ENDIF

	ELSE
	  PRINT *,'Unknown type  = ',
     $             INTRA_PROP,'!!!'
	  STOP
	ENDIF
	
        RETURN
	END


        FUNCTION MASS_ATOM(ATOM)
        CHARACTER*2 ATOM
        REAL*8 MASS_ATOM 


        IF(ATOM .EQ. 'C ') THEN
          MASS_ATOM = 12.010000
        ELSEIF(ATOM .EQ. 'H ') THEN
          MASS_ATOM = 1.008000
        ELSEIF(ATOM .EQ. 'O ')THEN
          MASS_ATOM = 15.999400
        ELSEIF(ATOM .EQ. 'OW') THEN
          MASS_ATOM = 15.999400
        ELSEIF(ATOM .EQ. 'HW') THEN
          MASS_ATOM =  1.008000
        ELSEIF(ATOM .EQ. 'N ') THEN
          MASS_ATOM = 14.010000
        ELSEIF(ATOM .EQ. 'Cl') THEN
          MASS_ATOM = 35.0
        ELSEIF(ATOM .EQ. 'Es') THEN
          MASS_ATOM = 0.0
        ELSEIF(ATOM .EQ. 'AR') THEN
          MASS_ATOM = 40.0
        ELSEIF(ATOM .EQ. 'Br') THEN
          MASS_ATOM = 79.0
        ELSEIF(ATOM .EQ. 'P ') THEN
          MASS_ATOM = 30.97376
        ELSEIF(ATOM .EQ. 'C3') THEN
          MASS_ATOM = 12.010000 + 3.0*1.008000
        ELSE
          print *,'unknown mass for atom : ',ATOM
          stop
        ENDIF 

        END

      subroutine cleanstring(filei)
      character*8 filei,fileo
      character*1 recname(8)
      logical ex(8)
      equivalence (recname(1),fileo)
c
      fileo=filei
c
c  clean from empty spaces
c
       do kk=1,8
         ex(kk)=.true.
         if (recname(kk).eq.' ') ex(kk)=.false.
       end do
c
       ns=0
       do kk=1,8
         if (ex(kk)) then
            ns=ns+1
            recname(ns)=recname(kk)
         end if
        end do
c
       do kk=ns+1,8
          recname(kk)=' '
       end do

c
       filei=fileo
       return
       end


