+!---------reading nucleic acid parameters for rotamers-------------------
+ allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
+ do i=1,ntyp_molec(2)
+ read (irotam_nucl,*)
+ do j=1,9
+ read(irotam_nucl,*) sc_parmin_nucl(j,i)
+ enddo
+ enddo
+ close(irotam_nucl)
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) "Base rotamer parameters"
+ do i=1,ntyp_molec(2)
+ write (iout,'(a)') restyp(i,2)
+ write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
+ enddo
+ endif
+
+
+ read (ifourier,*) nloctyp
+!el write(iout,*)"nloctyp",nloctyp
+ SPLIT_FOURIERTOR = nloctyp.lt.0
+ nloctyp = iabs(nloctyp)
+#ifdef NEWCORR
+ if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
+ print *,"shape",shape(itype2loc)
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(bnew1(3,2,-nloctyp:nloctyp))
+ allocate(bnew2(3,2,-nloctyp:nloctyp))
+ allocate(ccnew(3,2,-nloctyp:nloctyp))
+ allocate(ddnew(3,2,-nloctyp:nloctyp))
+ allocate(e0new(3,-nloctyp:nloctyp))
+ allocate(eenew(2,2,2,-nloctyp:nloctyp))
+ allocate(bnew1tor(3,2,-nloctyp:nloctyp))
+ allocate(bnew2tor(3,2,-nloctyp:nloctyp))
+ allocate(ccnewtor(3,2,-nloctyp:nloctyp))
+ allocate(ddnewtor(3,2,-nloctyp:nloctyp))
+ allocate(e0newtor(3,-nloctyp:nloctyp))
+ allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
+
+ read (ifourier,*) (itype2loc(i),i=1,ntyp)
+ read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
+ itype2loc(ntyp1)=nloctyp
+ iloctyp(nloctyp)=ntyp1
+ do i=1,ntyp1
+ itype2loc(-i)=-itype2loc(i)
+ enddo
+#else
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(itype2loc(-ntyp1:ntyp1))
+ iloctyp(0)=10
+ iloctyp(1)=9
+ iloctyp(2)=20
+ iloctyp(3)=ntyp1
+#endif
+ do i=1,nloctyp
+ iloctyp(-i)=-iloctyp(i)
+ enddo
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+#ifdef NEWCORR
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR",i
+ read (ifourier,*)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew1(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1"
+!c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew2(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2"
+!c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ccnew(kk,1,i)
+ read (ifourier,*) ccnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ddnew(kk,1,i)
+ read (ifourier,*) ddnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW"
+!c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*) eenew(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1"
+!c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*) e0new(ii,i)
+ enddo
+!c write (iout,*) (e0new(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnew(ii,1,i)=ccnew(ii,1,i)/2
+ ccnew(ii,2,i)=ccnew(ii,2,i)/2
+ ddnew(ii,1,i)=ddnew(ii,1,i)/2
+ ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1(ii,1,-i)= bnew1(ii,1,i)
+ bnew1(ii,2,-i)=-bnew1(ii,2,i)
+ bnew2(ii,1,-i)= bnew2(ii,1,i)
+ bnew2(ii,2,-i)=-bnew2(ii,2,i)
+ enddo
+ do ii=1,3
+!c ccnew(ii,1,i)=ccnew(ii,1,i)/2
+!c ccnew(ii,2,i)=ccnew(ii,2,i)/2
+!c ddnew(ii,1,i)=ddnew(ii,1,i)/2
+!c ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ ccnew(ii,1,-i)=ccnew(ii,1,i)
+ ccnew(ii,2,-i)=-ccnew(ii,2,i)
+ ddnew(ii,1,-i)=ddnew(ii,1,i)
+ ddnew(ii,2,-i)=-ddnew(ii,2,i)
+ enddo
+ e0new(1,-i)= e0new(1,i)
+ e0new(2,-i)=-e0new(2,i)
+ e0new(3,-i)=-e0new(3,i)
+ do kk=1,2
+ eenew(kk,1,1,-i)= eenew(kk,1,1,i)
+ eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
+ eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
+ eenew(kk,2,2,-i)= eenew(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') "Coefficients of the multibody terms"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of C"
+ write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of D"
+ write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of E"
+ write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ IF (SPLIT_FOURIERTOR) THEN
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR TOR",i
+ read (ifourier,*)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew1tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1 TOR"
+!c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew2tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2 TOR"
+!c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ccnewtor(kk,1,i)
+ read (ifourier,*) ccnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW TOR"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ddnewtor(kk,1,i)
+ read (ifourier,*) ddnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW TOR"
+!c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*) eenewtor(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1 TOR"
+!c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*) e0newtor(ii,i)
+ enddo
+!c write (iout,*) (e0newtor(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW TOR"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
+ ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
+ ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
+ ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
+ bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
+ bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
+ bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
+ ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
+ ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
+ ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
+ enddo
+ e0newtor(1,-i)= e0newtor(1,i)
+ e0newtor(2,-i)=-e0newtor(2,i)
+ e0newtor(3,-i)=-e0newtor(3,i)
+ do kk=1,2
+ eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
+ eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
+ eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
+ eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') &
+ "Single-body coefficients of the torsional potentials"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1tor"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') &
+ j,(bnew1tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2tor"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') &
+ j,(bnew2tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of Ctor"
+ write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Dtor"
+ write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Etor"
+ write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') &
+ j,k,(eenewtor(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ ELSE
+ do i=-nloctyp+1,nloctyp-1
+ do ii=1,3
+ do j=1,2
+ bnew1tor(ii,j,i)=bnew1(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ do j=1,2
+ bnew2tor(ii,j,i)=bnew2(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnew(ii,1,i)
+ ccnewtor(ii,2,i)=ccnew(ii,2,i)
+ ddnewtor(ii,1,i)=ddnew(ii,1,i)
+ ddnewtor(ii,2,i)=ddnew(ii,2,i)
+ enddo
+ enddo
+ ENDIF !SPLIT_FOURIER_TOR
+#else
+ allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(b(13,-nloctyp-1:nloctyp+1))
+ if (lprint) &
+ write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
+ do i=0,nloctyp-1
+ read (ifourier,*)
+ read (ifourier,*) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ if (i.gt.0) then
+ b(2,-i)= b(2,i)
+ b(3,-i)= b(3,i)
+ b(4,-i)=-b(4,i)
+ b(5,-i)=-b(5,i)
+ endif
+ CCold(1,1,i)= b(7,i)
+ CCold(2,2,i)=-b(7,i)
+ CCold(2,1,i)= b(9,i)
+ CCold(1,2,i)= b(9,i)
+ CCold(1,1,-i)= b(7,i)
+ CCold(2,2,-i)=-b(7,i)
+ CCold(2,1,-i)=-b(9,i)
+ CCold(1,2,-i)=-b(9,i)
+ DDold(1,1,i)= b(6,i)
+ DDold(2,2,i)=-b(6,i)
+ DDold(2,1,i)= b(8,i)
+ DDold(1,2,i)= b(8,i)
+ DDold(1,1,-i)= b(6,i)
+ DDold(2,2,-i)=-b(6,i)
+ DDold(2,1,-i)=-b(8,i)
+ DDold(1,2,-i)=-b(8,i)
+ EEold(1,1,i)= b(10,i)+b(11,i)
+ EEold(2,2,i)=-b(10,i)+b(11,i)
+ EEold(2,1,i)= b(12,i)-b(13,i)
+ EEold(1,2,i)= b(12,i)+b(13,i)
+ EEold(1,1,-i)= b(10,i)+b(11,i)
+ EEold(2,2,-i)=-b(10,i)+b(11,i)
+ EEold(2,1,-i)=-b(12,i)+b(13,i)
+ EEold(1,2,-i)=-b(12,i)-b(13,i)
+ write(iout,*) "TU DOCHODZE"
+ print *,"JESTEM"
+ enddo
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) &
+ "Coefficients of the cumulants (independent of valence angles)"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,*) 'B1'
+ write(iout,'(2f10.5)') B(3,i),B(5,i)
+ write (iout,*) 'B2'
+ write(iout,'(2f10.5)') B(2,i),B(4,i)
+ write (iout,*) 'CC'
+ do j=1,2
+ write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
+ enddo
+ write(iout,*) 'DD'
+ do j=1,2
+ write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
+ enddo
+ write(iout,*) 'EE'
+ do j=1,2
+ write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
+ enddo
+ enddo
+ endif
+#endif