8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 sigmasctube,krad2,ract
749 integer :: ichir1,ichir2,ijunk,irdiff
751 character*80 temp1,mychar
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
817 allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
818 allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) &
874 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
876 if (oldion.eq.1) then
878 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
879 print *,msc(i,5),restok(i,5)
883 read (iion,*) ncatprotparm
884 allocate(catprm(ncatprotparm,4))
886 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
890 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
894 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
897 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
898 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
902 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
903 restyp(i,5),"-",restyp(j,2)
906 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
907 !----------------------------------------------------
908 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
909 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
910 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
911 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
912 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
913 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
923 allocate(liptranene(ntyp))
924 !C reading lipid parameters
925 write (iout,*) "iliptranpar",iliptranpar
927 read(iliptranpar,*) pepliptran
930 read(iliptranpar,*) liptranene(i)
931 print *,liptranene(i)
934 ! water parmaters entalphy
935 allocate(awaterenta(0:400))
936 allocate(bwaterenta(0:400))
937 allocate(cwaterenta(0:400))
938 allocate(dwaterenta(0:400))
939 allocate(awaterentro(0:400))
940 allocate(bwaterentro(0:400))
941 allocate(cwaterentro(0:400))
942 allocate(dwaterentro(0:400))
944 read(iwaterwater,*) mychar
945 read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
946 cwaterenta(0),dwaterenta(0)
948 read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
949 cwaterenta(i),dwaterenta(i)
950 irdiff=int((ract-2.06d0)*50.0d0)+1
951 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
953 ! water parmaters entrophy
954 read(iwaterwater,*) mychar
955 read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
956 cwaterentro(0),dwaterentro(0)
958 read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
959 cwaterentro(i),dwaterentro(i)
960 irdiff=int((ract-2.06d0)*50.0d0)+1
961 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
967 ! Read the parameters of the probability distribution/energy expression
968 ! of the virtual-bond valence angles theta
971 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
972 (bthet(j,i,1,1),j=1,2)
973 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
974 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
975 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
979 athet(1,i,1,-1)=athet(1,i,1,1)
980 athet(2,i,1,-1)=athet(2,i,1,1)
981 bthet(1,i,1,-1)=-bthet(1,i,1,1)
982 bthet(2,i,1,-1)=-bthet(2,i,1,1)
983 athet(1,i,-1,1)=-athet(1,i,1,1)
984 athet(2,i,-1,1)=-athet(2,i,1,1)
985 bthet(1,i,-1,1)=bthet(1,i,1,1)
986 bthet(2,i,-1,1)=bthet(2,i,1,1)
990 athet(1,i,-1,-1)=athet(1,-i,1,1)
991 athet(2,i,-1,-1)=-athet(2,-i,1,1)
992 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
993 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
994 athet(1,i,-1,1)=athet(1,-i,1,1)
995 athet(2,i,-1,1)=-athet(2,-i,1,1)
996 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
997 bthet(2,i,-1,1)=bthet(2,-i,1,1)
998 athet(1,i,1,-1)=-athet(1,-i,1,1)
999 athet(2,i,1,-1)=athet(2,-i,1,1)
1000 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1001 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1002 theta0(i)=theta0(-i)
1006 polthet(j,i)=polthet(j,-i)
1009 gthet(j,i)=gthet(j,-i)
1015 if (.not.LaTeX) then
1016 write (iout,'(a)') &
1017 'Parameters of the virtual-bond valence angles:'
1018 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
1019 ' ATHETA0 ',' A1 ',' A2 ',&
1022 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1023 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
1025 write (iout,'(/a/9x,5a/79(1h-))') &
1026 'Parameters of the expression for sigma(theta_c):',&
1027 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
1028 ' ALPH3 ',' SIGMA0C '
1030 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1031 (polthet(j,i),j=0,3),sigc0(i)
1033 write (iout,'(/a/9x,5a/79(1h-))') &
1034 'Parameters of the second gaussian:',&
1035 ' THETA0 ',' SIGMA0 ',' G1 ',&
1038 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1039 sig0(i),(gthet(j,i),j=1,3)
1042 write (iout,'(a)') &
1043 'Parameters of the virtual-bond valence angles:'
1044 write (iout,'(/a/9x,5a/79(1h-))') &
1045 'Coefficients of expansion',&
1046 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1047 ' b1*10^1 ',' b2*10^1 '
1049 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1050 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1051 (10*bthet(j,i,1,1),j=1,2)
1053 write (iout,'(/a/9x,5a/79(1h-))') &
1054 'Parameters of the expression for sigma(theta_c):',&
1055 ' alpha0 ',' alph1 ',' alph2 ',&
1056 ' alhp3 ',' sigma0c '
1058 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1059 (polthet(j,i),j=0,3),sigc0(i)
1061 write (iout,'(/a/9x,5a/79(1h-))') &
1062 'Parameters of the second gaussian:',&
1063 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1066 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1067 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1073 ! Read the parameters of Utheta determined from ab initio surfaces
1074 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1076 IF (tor_mode.eq.0) THEN
1077 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1078 ntheterm3,nsingle,ndouble
1079 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1081 !----------------------------------------------------
1082 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1083 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1084 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1085 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1086 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1087 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1088 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1089 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1090 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1091 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1092 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1093 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1094 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1095 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1096 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1097 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1098 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1099 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1100 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1101 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1102 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1103 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1104 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1105 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1107 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1109 ithetyp(i)=-ithetyp(-i)
1112 aa0thet(:,:,:,:)=0.0d0
1113 aathet(:,:,:,:,:)=0.0d0
1114 bbthet(:,:,:,:,:,:)=0.0d0
1115 ccthet(:,:,:,:,:,:)=0.0d0
1116 ddthet(:,:,:,:,:,:)=0.0d0
1117 eethet(:,:,:,:,:,:)=0.0d0
1118 ffthet(:,:,:,:,:,:,:)=0.0d0
1119 ggthet(:,:,:,:,:,:,:)=0.0d0
1121 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1123 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1124 ! VAR:1=non-glicyne non-proline 2=proline
1125 ! VAR:negative values for D-aminoacid
1127 do j=-nthetyp,nthetyp
1128 do k=-nthetyp,nthetyp
1129 read (ithep,'(6a)',end=111,err=111) res1
1130 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1131 ! VAR: aa0thet is variable describing the average value of Foureir
1132 ! VAR: expansion series
1133 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1134 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1135 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1136 read (ithep,*,end=111,err=111) &
1137 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1138 read (ithep,*,end=111,err=111) &
1139 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1140 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1141 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1142 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1144 read (ithep,*,end=111,err=111) &
1145 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1146 ffthet(lll,llll,ll,i,j,k,iblock),&
1147 ggthet(llll,lll,ll,i,j,k,iblock),&
1148 ggthet(lll,llll,ll,i,j,k,iblock),&
1149 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1154 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1155 ! coefficients of theta-and-gamma-dependent terms are zero.
1156 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1157 ! RECOMENTDED AFTER VERSION 3.3)
1161 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1162 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1164 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1165 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1168 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1170 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1173 ! AND COMMENT THE LOOPS BELOW
1177 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1178 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1180 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1181 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1184 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1186 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1191 ! Substitution for D aminoacids from symmetry.
1194 do j=-nthetyp,nthetyp
1195 do k=-nthetyp,nthetyp
1196 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1198 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1202 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1203 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1204 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1205 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1211 ffthet(llll,lll,ll,i,j,k,iblock)= &
1212 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1213 ffthet(lll,llll,ll,i,j,k,iblock)= &
1214 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1215 ggthet(llll,lll,ll,i,j,k,iblock)= &
1216 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1217 ggthet(lll,llll,ll,i,j,k,iblock)= &
1218 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1227 ! Control printout of the coefficients of virtual-bond-angle potentials
1230 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1235 write (iout,'(//4a)') &
1236 'Type ',onelett(i),onelett(j),onelett(k)
1237 write (iout,'(//a,10x,a)') " l","a[l]"
1238 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1239 write (iout,'(i2,1pe15.5)') &
1240 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1242 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1243 "b",l,"c",l,"d",l,"e",l
1245 write (iout,'(i2,4(1pe15.5))') m,&
1246 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1247 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1251 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1252 "f+",l,"f-",l,"g+",l,"g-",l
1255 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1256 ffthet(n,m,l,i,j,k,iblock),&
1257 ffthet(m,n,l,i,j,k,iblock),&
1258 ggthet(n,m,l,i,j,k,iblock),&
1259 ggthet(m,n,l,i,j,k,iblock)
1270 !C here will be the apropriate recalibrating for D-aminoacid
1271 read (ithep,*,end=121,err=121) nthetyp
1272 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1273 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1274 do i=-nthetyp+1,nthetyp-1
1275 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1276 do j=0,nbend_kcc_Tb(i)
1277 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1281 write (iout,'(a)') &
1282 "Parameters of the valence-only potentials"
1283 do i=-nthetyp+1,nthetyp-1
1284 write (iout,'(2a)') "Type ",toronelet(i)
1285 do k=0,nbend_kcc_Tb(i)
1286 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1292 write (2,*) "Start reading THETA_PDB",ithep_pdb
1294 ! write (2,*) 'i=',i
1295 read (ithep_pdb,*,err=111,end=111) &
1296 a0thet(i),(athet(j,i,1,1),j=1,2),&
1297 (bthet(j,i,1,1),j=1,2)
1298 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1299 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1300 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1301 sigc0(i)=sigc0(i)**2
1304 athet(1,i,1,-1)=athet(1,i,1,1)
1305 athet(2,i,1,-1)=athet(2,i,1,1)
1306 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1307 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1308 athet(1,i,-1,1)=-athet(1,i,1,1)
1309 athet(2,i,-1,1)=-athet(2,i,1,1)
1310 bthet(1,i,-1,1)=bthet(1,i,1,1)
1311 bthet(2,i,-1,1)=bthet(2,i,1,1)
1314 a0thet(i)=a0thet(-i)
1315 athet(1,i,-1,-1)=athet(1,-i,1,1)
1316 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1317 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1318 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1319 athet(1,i,-1,1)=athet(1,-i,1,1)
1320 athet(2,i,-1,1)=-athet(2,-i,1,1)
1321 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1322 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1323 athet(1,i,1,-1)=-athet(1,-i,1,1)
1324 athet(2,i,1,-1)=athet(2,-i,1,1)
1325 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1326 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1327 theta0(i)=theta0(-i)
1331 polthet(j,i)=polthet(j,-i)
1334 gthet(j,i)=gthet(j,-i)
1337 write (2,*) "End reading THETA_PDB"
1341 !--------------- Reading theta parameters for nucleic acid-------
1342 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1343 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1344 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1345 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1346 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1347 nthetyp_nucl+1,nthetyp_nucl+1))
1348 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1349 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1350 nthetyp_nucl+1,nthetyp_nucl+1))
1351 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1352 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1353 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1354 nthetyp_nucl+1,nthetyp_nucl+1))
1355 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1356 nthetyp_nucl+1,nthetyp_nucl+1))
1357 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1358 nthetyp_nucl+1,nthetyp_nucl+1))
1359 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1360 nthetyp_nucl+1,nthetyp_nucl+1))
1361 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1362 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1363 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1364 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1365 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1366 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1368 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1369 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1371 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1373 aa0thet_nucl(:,:,:)=0.0d0
1374 aathet_nucl(:,:,:,:)=0.0d0
1375 bbthet_nucl(:,:,:,:,:)=0.0d0
1376 ccthet_nucl(:,:,:,:,:)=0.0d0
1377 ddthet_nucl(:,:,:,:,:)=0.0d0
1378 eethet_nucl(:,:,:,:,:)=0.0d0
1379 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1380 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1385 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1386 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1387 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1388 read (ithep_nucl,*,end=111,err=111) &
1389 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1390 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1391 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1392 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1393 read (ithep_nucl,*,end=111,err=111) &
1394 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1395 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1396 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1401 !-------------------------------------------
1402 allocate(nlob(ntyp1)) !(ntyp1)
1403 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1404 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1405 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1412 gaussc(:,:,:,:)=0.0D0
1416 ! Read the parameters of the probability distribution/energy expression
1417 ! of the side chains.
1420 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1424 dsc_inv(i)=1.0D0/dsc(i)
1435 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1436 ((blower(k,l,1),l=1,k),k=1,3)
1437 censc(1,1,-i)=censc(1,1,i)
1438 censc(2,1,-i)=censc(2,1,i)
1439 censc(3,1,-i)=-censc(3,1,i)
1441 read (irotam,*,end=112,err=112) bsc(j,i)
1442 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1443 ((blower(k,l,j),l=1,k),k=1,3)
1444 censc(1,j,-i)=censc(1,j,i)
1445 censc(2,j,-i)=censc(2,j,i)
1446 censc(3,j,-i)=-censc(3,j,i)
1447 ! BSC is amplitude of Gaussian
1454 akl=akl+blower(k,m,j)*blower(l,m,j)
1458 if (((k.eq.3).and.(l.ne.3)) &
1459 .or.((l.eq.3).and.(k.ne.3))) then
1460 gaussc(k,l,j,-i)=-akl
1461 gaussc(l,k,j,-i)=-akl
1463 gaussc(k,l,j,-i)=akl
1464 gaussc(l,k,j,-i)=akl
1473 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1476 if (nlobi.gt.0) then
1478 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1479 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1480 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1481 'log h',(bsc(j,i),j=1,nlobi)
1482 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1483 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1485 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1486 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1489 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1490 write (iout,'(a,f10.4,4(16x,f10.4))') &
1491 'Center ',(bsc(j,i),j=1,nlobi)
1492 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1501 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1502 ! added by Urszula Kozlowska 07/11/2007
1504 !el Maximum number of SC local term fitting function coefficiants
1505 !el integer,parameter :: maxsccoef=65
1507 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1510 read (irotam,*,end=112,err=112)
1512 read (irotam,*,end=112,err=112)
1515 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1519 !---------reading nucleic acid parameters for rotamers-------------------
1520 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1521 do i=1,ntyp_molec(2)
1522 read (irotam_nucl,*,end=112,err=112)
1524 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1530 write (iout,*) "Base rotamer parameters"
1531 do i=1,ntyp_molec(2)
1532 write (iout,'(a)') restyp(i,2)
1533 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1538 ! Read the parameters of the probability distribution/energy expression
1539 ! of the side chains.
1541 write (2,*) "Start reading ROTAM_PDB"
1543 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1547 dsc_inv(i)=1.0D0/dsc(i)
1558 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1559 ((blower(k,l,1),l=1,k),k=1,3)
1561 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1562 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1563 ((blower(k,l,j),l=1,k),k=1,3)
1570 akl=akl+blower(k,m,j)*blower(l,m,j)
1580 write (2,*) "End reading ROTAM_PDB"
1586 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1587 !C interaction energy of the Gly, Ala, and Pro prototypes.
1589 read (ifourier,*) nloctyp
1590 SPLIT_FOURIERTOR = nloctyp.lt.0
1591 nloctyp = iabs(nloctyp)
1592 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1593 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1594 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1595 !C allocate(ctilde(2,2,nres))
1596 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1597 !C allocate(gtb1(2,nres))
1598 !C allocate(gtb2(2,nres))
1599 !C allocate(cc(2,2,nres))
1600 !C allocate(dd(2,2,nres))
1601 !C allocate(ee(2,2,nres))
1602 !C allocate(gtcc(2,2,nres))
1603 !C allocate(gtdd(2,2,nres))
1604 !C allocate(gtee(2,2,nres))
1607 allocate(itype2loc(-ntyp1:ntyp1))
1608 allocate(iloctyp(-nloctyp:nloctyp))
1609 allocate(bnew1(3,2,-nloctyp:nloctyp))
1610 allocate(bnew2(3,2,-nloctyp:nloctyp))
1611 allocate(ccnew(3,2,-nloctyp:nloctyp))
1612 allocate(ddnew(3,2,-nloctyp:nloctyp))
1613 allocate(e0new(3,-nloctyp:nloctyp))
1614 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1615 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1616 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1617 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1618 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1619 allocate(e0newtor(3,-nloctyp:nloctyp))
1620 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1633 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1634 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1635 itype2loc(ntyp1)=nloctyp
1636 iloctyp(nloctyp)=ntyp1
1638 itype2loc(-i)=-itype2loc(i)
1641 allocate(iloctyp(-nloctyp:nloctyp))
1642 allocate(itype2loc(-ntyp1:ntyp1))
1649 iloctyp(-i)=-iloctyp(i)
1651 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1652 !c write (iout,*) "nloctyp",nloctyp,
1653 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1654 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1655 !c write (iout,*) "nloctyp",nloctyp,
1656 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1659 !c write (iout,*) "NEWCORR",i
1660 read (ifourier,*,end=115,err=115)
1663 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1666 !c write (iout,*) "NEWCORR BNEW1"
1667 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1670 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1673 !c write (iout,*) "NEWCORR BNEW2"
1674 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1676 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1677 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1679 !c write (iout,*) "NEWCORR CCNEW"
1680 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1682 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1683 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1685 !c write (iout,*) "NEWCORR DDNEW"
1686 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1690 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1694 !c write (iout,*) "NEWCORR EENEW1"
1695 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1697 read (ifourier,*,end=115,err=115) e0new(ii,i)
1699 !c write (iout,*) (e0new(ii,i),ii=1,3)
1701 !c write (iout,*) "NEWCORR EENEW"
1704 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1705 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1706 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1707 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1712 bnew1(ii,1,-i)= bnew1(ii,1,i)
1713 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1714 bnew2(ii,1,-i)= bnew2(ii,1,i)
1715 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1718 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1719 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1720 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1721 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1722 ccnew(ii,1,-i)=ccnew(ii,1,i)
1723 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1724 ddnew(ii,1,-i)=ddnew(ii,1,i)
1725 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1727 e0new(1,-i)= e0new(1,i)
1728 e0new(2,-i)=-e0new(2,i)
1729 e0new(3,-i)=-e0new(3,i)
1731 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1732 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1733 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1734 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1738 write (iout,'(a)') "Coefficients of the multibody terms"
1739 do i=-nloctyp+1,nloctyp-1
1740 write (iout,*) "Type: ",onelet(iloctyp(i))
1741 write (iout,*) "Coefficients of the expansion of B1"
1743 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1745 write (iout,*) "Coefficients of the expansion of B2"
1747 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1749 write (iout,*) "Coefficients of the expansion of C"
1750 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1751 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1752 write (iout,*) "Coefficients of the expansion of D"
1753 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1754 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1755 write (iout,*) "Coefficients of the expansion of E"
1756 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1759 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1764 IF (SPLIT_FOURIERTOR) THEN
1766 !c write (iout,*) "NEWCORR TOR",i
1767 read (ifourier,*,end=115,err=115)
1770 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1773 !c write (iout,*) "NEWCORR BNEW1 TOR"
1774 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1777 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1780 !c write (iout,*) "NEWCORR BNEW2 TOR"
1781 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1783 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1784 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1786 !c write (iout,*) "NEWCORR CCNEW TOR"
1787 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1789 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1790 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1792 !c write (iout,*) "NEWCORR DDNEW TOR"
1793 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1797 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1801 !c write (iout,*) "NEWCORR EENEW1 TOR"
1802 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1804 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1806 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1808 !c write (iout,*) "NEWCORR EENEW TOR"
1811 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1812 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1813 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1814 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1819 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1820 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1821 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1822 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1825 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1826 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1827 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1828 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1830 e0newtor(1,-i)= e0newtor(1,i)
1831 e0newtor(2,-i)=-e0newtor(2,i)
1832 e0newtor(3,-i)=-e0newtor(3,i)
1834 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1835 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1836 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1837 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1841 write (iout,'(a)') &
1842 "Single-body coefficients of the torsional potentials"
1843 do i=-nloctyp+1,nloctyp-1
1844 write (iout,*) "Type: ",onelet(iloctyp(i))
1845 write (iout,*) "Coefficients of the expansion of B1tor"
1847 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1848 j,(bnew1tor(k,j,i),k=1,3)
1850 write (iout,*) "Coefficients of the expansion of B2tor"
1852 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1853 j,(bnew2tor(k,j,i),k=1,3)
1855 write (iout,*) "Coefficients of the expansion of Ctor"
1856 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1857 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1858 write (iout,*) "Coefficients of the expansion of Dtor"
1859 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1860 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1861 write (iout,*) "Coefficients of the expansion of Etor"
1862 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1865 write (iout,'(1hE,2i1,2f10.5)') &
1866 j,k,(eenewtor(l,j,k,i),l=1,2)
1872 do i=-nloctyp+1,nloctyp-1
1875 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1880 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1884 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1885 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1886 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1887 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1890 ENDIF !SPLIT_FOURIER_TOR
1892 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1893 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1894 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1895 allocate(b(13,-nloctyp-1:nloctyp+1))
1897 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1899 read (ifourier,*,end=115,err=115)
1900 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1902 write (iout,*) 'Type ',onelet(iloctyp(i))
1903 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1917 !c B1tilde(1,i) = b(3)
1918 !c! B1tilde(2,i) =-b(5)
1919 !c! B1tilde(1,-i) =-b(3)
1920 !c! B1tilde(2,-i) =b(5)
1921 !c! b1tilde(1,i)=0.0d0
1922 !c b1tilde(2,i)=0.0d0
1927 !cc B1tilde(1,i) = b(3,i)
1928 !cc B1tilde(2,i) =-b(5,i)
1929 !c B1tilde(1,-i) =-b(3,i)
1930 !c B1tilde(2,-i) =b(5,i)
1931 !cc b1tilde(1,i)=0.0d0
1932 !cc b1tilde(2,i)=0.0d0
1933 !cc B2(1,i) = b(2,i)
1934 !cc B2(2,i) = b(4,i)
1936 !c B2(2,-i) =-b(4,i)
1940 CCold(1,1,i)= b(7,i)
1941 CCold(2,2,i)=-b(7,i)
1942 CCold(2,1,i)= b(9,i)
1943 CCold(1,2,i)= b(9,i)
1944 CCold(1,1,-i)= b(7,i)
1945 CCold(2,2,-i)=-b(7,i)
1946 CCold(2,1,-i)=-b(9,i)
1947 CCold(1,2,-i)=-b(9,i)
1952 !c Ctilde(1,1,i)= CCold(1,1,i)
1953 !c Ctilde(1,2,i)= CCold(1,2,i)
1954 !c Ctilde(2,1,i)=-CCold(2,1,i)
1955 !c Ctilde(2,2,i)=-CCold(2,2,i)
1960 !c Ctilde(1,1,i)= CCold(1,1,i)
1961 !c Ctilde(1,2,i)= CCold(1,2,i)
1962 !c Ctilde(2,1,i)=-CCold(2,1,i)
1963 !c Ctilde(2,2,i)=-CCold(2,2,i)
1965 !c Ctilde(1,1,i)=0.0d0
1966 !c Ctilde(1,2,i)=0.0d0
1967 !c Ctilde(2,1,i)=0.0d0
1968 !c Ctilde(2,2,i)=0.0d0
1969 DDold(1,1,i)= b(6,i)
1970 DDold(2,2,i)=-b(6,i)
1971 DDold(2,1,i)= b(8,i)
1972 DDold(1,2,i)= b(8,i)
1973 DDold(1,1,-i)= b(6,i)
1974 DDold(2,2,-i)=-b(6,i)
1975 DDold(2,1,-i)=-b(8,i)
1976 DDold(1,2,-i)=-b(8,i)
1981 !c Dtilde(1,1,i)= DD(1,1,i)
1982 !c Dtilde(1,2,i)= DD(1,2,i)
1983 !c Dtilde(2,1,i)=-DD(2,1,i)
1984 !c Dtilde(2,2,i)=-DD(2,2,i)
1986 !c Dtilde(1,1,i)=0.0d0
1987 !c Dtilde(1,2,i)=0.0d0
1988 !c Dtilde(2,1,i)=0.0d0
1989 !c Dtilde(2,2,i)=0.0d0
1990 EEold(1,1,i)= b(10,i)+b(11,i)
1991 EEold(2,2,i)=-b(10,i)+b(11,i)
1992 EEold(2,1,i)= b(12,i)-b(13,i)
1993 EEold(1,2,i)= b(12,i)+b(13,i)
1994 EEold(1,1,-i)= b(10,i)+b(11,i)
1995 EEold(2,2,-i)=-b(10,i)+b(11,i)
1996 EEold(2,1,-i)=-b(12,i)+b(13,i)
1997 EEold(1,2,-i)=-b(12,i)-b(13,i)
1998 write(iout,*) "TU DOCHODZE"
2004 !c ee(2,1,i)=ee(1,2,i)
2009 "Coefficients of the cumulants (independent of valence angles)"
2010 do i=-nloctyp+1,nloctyp-1
2011 write (iout,*) 'Type ',onelet(iloctyp(i))
2013 write(iout,'(2f10.5)') B(3,i),B(5,i)
2015 write(iout,'(2f10.5)') B(2,i),B(4,i)
2018 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
2022 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
2026 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
2035 ! Read torsional parameters in old format
2037 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
2039 read (itorp,*,end=113,err=113) ntortyp,nterm_old
2040 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
2041 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2043 !el from energy module--------
2044 allocate(v1(nterm_old,ntortyp,ntortyp))
2045 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2046 !el---------------------------
2051 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2057 write (iout,'(/a/)') 'Torsional constants:'
2060 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2061 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2067 ! Read torsional parameters
2069 IF (TOR_MODE.eq.0) THEN
2070 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2071 read (itorp,*,end=113,err=113) ntortyp
2072 !el from energy module---------
2073 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2074 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2076 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2077 allocate(vlor2(maxlor,ntortyp,ntortyp))
2078 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2079 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2081 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2082 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2083 !el---------------------------
2086 !el---------------------------
2088 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2090 itortyp(i)=-itortyp(-i)
2092 itortyp(ntyp1)=ntortyp
2093 itortyp(-ntyp1)=-ntortyp
2095 write (iout,*) 'ntortyp',ntortyp
2097 do j=-ntortyp+1,ntortyp-1
2098 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2100 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2101 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2104 do k=1,nterm(i,j,iblock)
2105 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2107 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2108 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2109 v0ij=v0ij+si*v1(k,i,j,iblock)
2111 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2112 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2113 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2115 do k=1,nlor(i,j,iblock)
2116 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2117 vlor2(k,i,j),vlor3(k,i,j)
2118 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2121 v0(-i,-j,iblock)=v0ij
2127 write (iout,'(/a/)') 'Torsional constants:'
2129 do i=-ntortyp,ntortyp
2130 do j=-ntortyp,ntortyp
2131 write (iout,*) 'ityp',i,' jtyp',j
2132 write (iout,*) 'Fourier constants'
2133 do k=1,nterm(i,j,iblock)
2134 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2137 write (iout,*) 'Lorenz constants'
2138 do k=1,nlor(i,j,iblock)
2139 write (iout,'(3(1pe15.5))') &
2140 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2146 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2148 ! 6/23/01 Read parameters for double torsionals
2150 !el from energy module------------
2151 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2152 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2153 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2154 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2155 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2156 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2157 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2158 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2159 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2160 !---------------------------------
2164 do j=-ntortyp+1,ntortyp-1
2165 do k=-ntortyp+1,ntortyp-1
2166 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2167 ! write (iout,*) "OK onelett",
2170 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2171 .or. t3.ne.toronelet(k)) then
2172 write (iout,*) "Error in double torsional parameter file",&
2175 call MPI_Finalize(Ierror)
2177 stop "Error in double torsional parameter file"
2179 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2180 ntermd_2(i,j,k,iblock)
2181 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2182 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2183 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2184 ntermd_1(i,j,k,iblock))
2185 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2186 ntermd_1(i,j,k,iblock))
2187 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2188 ntermd_1(i,j,k,iblock))
2189 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2190 ntermd_1(i,j,k,iblock))
2191 ! Martix of D parameters for one dimesional foureir series
2192 do l=1,ntermd_1(i,j,k,iblock)
2193 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2194 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2195 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2196 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2197 ! write(iout,*) "whcodze" ,
2198 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2200 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2201 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2202 v2s(m,l,i,j,k,iblock),&
2203 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2204 ! Martix of D parameters for two dimesional fourier series
2205 do l=1,ntermd_2(i,j,k,iblock)
2207 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2208 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2209 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2210 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2219 write (iout,*) 'Constants for double torsionals'
2222 do j=-ntortyp+1,ntortyp-1
2223 do k=-ntortyp+1,ntortyp-1
2224 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2225 ' nsingle',ntermd_1(i,j,k,iblock),&
2226 ' ndouble',ntermd_2(i,j,k,iblock)
2228 write (iout,*) 'Single angles:'
2229 do l=1,ntermd_1(i,j,k,iblock)
2230 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2231 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2232 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2233 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2236 write (iout,*) 'Pairs of angles:'
2237 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2238 do l=1,ntermd_2(i,j,k,iblock)
2239 write (iout,'(i5,20f10.5)') &
2240 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2243 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2244 do l=1,ntermd_2(i,j,k,iblock)
2245 write (iout,'(i5,20f10.5)') &
2246 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2247 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2257 itype2loc(i)=itortyp(i)
2261 ELSE IF (TOR_MODE.eq.1) THEN
2263 !C read valence-torsional parameters
2264 read (itorp,*,end=121,err=121) ntortyp
2266 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2268 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2269 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2272 itype2loc(i)=itortyp(i)
2276 itortyp(i)=-itortyp(-i)
2278 do i=-ntortyp+1,ntortyp-1
2279 do j=-ntortyp+1,ntortyp-1
2280 !C first we read the cos and sin gamma parameters
2281 read (itorp,'(13x,a)',end=121,err=121) string
2282 write (iout,*) i,j,string
2283 read (itorp,*,end=121,err=121) &
2284 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2285 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2286 do k=1,nterm_kcc(j,i)
2287 do l=1,nterm_kcc_Tb(j,i)
2288 do ll=1,nterm_kcc_Tb(j,i)
2289 read (itorp,*,end=121,err=121) ii,jj,kk, &
2290 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2298 !c AL 4/8/16: Calculate coefficients from one-body parameters
2300 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2301 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2302 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2303 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2304 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2307 print *,i,itortyp(i)
2308 itortyp(i)=itype2loc(i)
2311 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2313 do i=-ntortyp+1,ntortyp-1
2314 do j=-ntortyp+1,ntortyp-1
2317 do k=1,nterm_kcc_Tb(j,i)
2318 do l=1,nterm_kcc_Tb(j,i)
2319 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2320 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2321 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2322 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2325 do k=1,nterm_kcc_Tb(j,i)
2326 do l=1,nterm_kcc_Tb(j,i)
2328 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2329 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2330 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2331 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2333 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2334 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2335 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2336 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2340 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2344 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2349 if (tor_mode.gt.0 .and. lprint) then
2350 !c Print valence-torsional parameters
2351 write (iout,'(a)') &
2352 "Parameters of the valence-torsional potentials"
2353 do i=-ntortyp+1,ntortyp-1
2354 do j=-ntortyp+1,ntortyp-1
2355 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2356 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2357 do k=1,nterm_kcc(j,i)
2358 do l=1,nterm_kcc_Tb(j,i)
2359 do ll=1,nterm_kcc_Tb(j,i)
2360 write (iout,'(3i5,2f15.4)')&
2361 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2369 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2370 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2371 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2372 !el from energy module---------
2373 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2374 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2376 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2377 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2378 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2379 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2381 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2382 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2383 !el---------------------------
2386 !el--------------------
2387 read (itorp_nucl,*,end=113,err=113) &
2388 (itortyp_nucl(i),i=1,ntyp_molec(2))
2389 ! print *,itortyp_nucl(:)
2390 !c write (iout,*) 'ntortyp',ntortyp
2393 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2394 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2397 do k=1,nterm_nucl(i,j)
2398 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2399 v0ij=v0ij+si*v1_nucl(k,i,j)
2402 do k=1,nlor_nucl(i,j)
2403 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2404 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2405 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2411 ! Read of Side-chain backbone correlation parameters
2412 ! Modified 11 May 2012 by Adasko
2415 read (isccor,*,end=119,err=119) nsccortyp
2417 !el from module energy-------------
2418 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2419 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2420 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2421 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2422 !-----------------------------------
2424 !el from module energy-------------
2425 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2427 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2429 isccortyp(i)=-isccortyp(-i)
2431 iscprol=isccortyp(20)
2432 ! write (iout,*) 'ntortyp',ntortyp
2434 !c maxinter is maximum interaction sites
2435 !el from module energy---------
2436 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2437 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2438 -nsccortyp:nsccortyp))
2439 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2440 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2441 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2442 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2443 !-----------------------------------
2447 read (isccor,*,end=119,err=119) &
2448 nterm_sccor(i,j),nlor_sccor(i,j)
2454 nterm_sccor(-i,j)=nterm_sccor(i,j)
2455 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2456 nterm_sccor(i,-j)=nterm_sccor(i,j)
2457 do k=1,nterm_sccor(i,j)
2458 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2460 if (j.eq.iscprol) then
2461 if (i.eq.isccortyp(10)) then
2462 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2463 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2465 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2466 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2467 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2468 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2469 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2470 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2471 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2472 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2475 if (i.eq.isccortyp(10)) then
2476 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2477 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2479 if (j.eq.isccortyp(10)) then
2480 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2481 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2483 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2484 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2485 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2486 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2487 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2488 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2492 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2493 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2494 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2495 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2498 do k=1,nlor_sccor(i,j)
2499 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2500 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2501 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2502 (1+vlor3sccor(k,i,j)**2)
2504 v0sccor(l,i,j)=v0ijsccor
2505 v0sccor(l,-i,j)=v0ijsccor1
2506 v0sccor(l,i,-j)=v0ijsccor2
2507 v0sccor(l,-i,-j)=v0ijsccor3
2513 !el from module energy-------------
2514 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2516 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2517 ! write (iout,*) 'ntortyp',ntortyp
2519 !c maxinter is maximum interaction sites
2520 !el from module energy---------
2521 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2522 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2523 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2524 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2525 !-----------------------------------
2529 read (isccor,*,end=119,err=119) &
2530 nterm_sccor(i,j),nlor_sccor(i,j)
2534 do k=1,nterm_sccor(i,j)
2535 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2537 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2540 do k=1,nlor_sccor(i,j)
2541 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2542 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2543 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2544 (1+vlor3sccor(k,i,j)**2)
2546 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2555 write (iout,'(/a/)') 'Torsional constants:'
2558 write (iout,*) 'ityp',i,' jtyp',j
2559 write (iout,*) 'Fourier constants'
2560 do k=1,nterm_sccor(i,j)
2561 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2563 write (iout,*) 'Lorenz constants'
2564 do k=1,nlor_sccor(i,j)
2565 write (iout,'(3(1pe15.5))') &
2566 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2573 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2574 ! interaction energy of the Gly, Ala, and Pro prototypes.
2577 ! Read electrostatic-interaction parameters
2582 write (iout,'(/a)') 'Electrostatic interaction constants:'
2583 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2584 'IT','JT','APP','BPP','AEL6','AEL3'
2586 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2587 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2588 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2589 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2594 app (i,j)=epp(i,j)*rri*rri
2595 bpp (i,j)=-2.0D0*epp(i,j)*rri
2596 ael6(i,j)=elpp6(i,j)*4.2D0**6
2597 ael3(i,j)=elpp3(i,j)*4.2D0**3
2599 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2605 ! Read side-chain interaction parameters.
2607 !el from module energy - COMMON.INTERACT-------
2608 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2609 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2610 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2612 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2613 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2614 allocate(epslip(ntyp,ntyp))
2622 !--------------------------------
2624 read (isidep,*,end=117,err=117) ipot,expon
2625 if (ipot.lt.1 .or. ipot.gt.5) then
2626 write (iout,'(2a)') 'Error while reading SC interaction',&
2627 'potential file - unknown potential type.'
2629 call MPI_Finalize(Ierror)
2635 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2636 ', exponents are ',expon,2*expon
2637 ! goto (10,20,30,30,40) ipot
2639 !----------------------- LJ potential ---------------------------------
2641 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2642 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2643 (sigma0(i),i=1,ntyp)
2645 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2646 write (iout,'(a/)') 'The epsilon array:'
2647 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2648 write (iout,'(/a)') 'One-body parameters:'
2649 write (iout,'(a,4x,a)') 'residue','sigma'
2650 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2653 !----------------------- LJK potential --------------------------------
2655 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2656 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2657 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2659 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2660 write (iout,'(a/)') 'The epsilon array:'
2661 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2662 write (iout,'(/a)') 'One-body parameters:'
2663 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2664 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2668 !---------------------- GB or BP potential -----------------------------
2671 ! print *,"I AM in SCELE",scelemode
2672 if (scelemode.eq.0) then
2674 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2676 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2677 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2678 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2679 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2681 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2684 ! For the GB potential convert sigma'**2 into chi'
2687 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2691 write (iout,'(/a/)') 'Parameters of the BP potential:'
2692 write (iout,'(a/)') 'The epsilon array:'
2693 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2694 write (iout,'(/a)') 'One-body parameters:'
2695 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2697 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2698 chip(i),alp(i),i=1,ntyp)
2701 ! print *,ntyp,"NTYP"
2702 allocate(icharge(ntyp1))
2703 ! print *,ntyp,icharge(i)
2705 read (isidep,*) (icharge(i),i=1,ntyp)
2706 print *,ntyp,icharge(i)
2707 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2708 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2709 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2710 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2711 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2712 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2713 allocate(epsintab(ntyp,ntyp))
2714 allocate(dtail(2,ntyp,ntyp))
2715 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2716 allocate(wqdip(2,ntyp,ntyp))
2717 allocate(wstate(4,ntyp,ntyp))
2718 allocate(dhead(2,2,ntyp,ntyp))
2719 allocate(nstate(ntyp,ntyp))
2720 allocate(debaykap(ntyp,ntyp))
2722 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2723 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2727 ! write (*,*) "Im in ALAB", i, " ", j
2729 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2730 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2731 chis(i,j),chis(j,i), & !2 w tej linii
2732 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2733 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2734 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2735 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2736 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2737 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2738 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2739 IF ((LaTeX).and.(i.gt.24)) then
2740 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2741 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2742 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2743 chis(i,j),chis(j,i) !2 w tej linii
2745 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2749 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2753 IF ((LaTeX).and.(i.gt.24)) then
2754 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2755 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2756 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2757 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2758 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2759 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2760 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2767 sigma(i,j) = sigma(j,i)
2768 sigmap1(i,j)=sigmap1(j,i)
2769 sigmap2(i,j)=sigmap2(j,i)
2770 sigiso1(i,j)=sigiso1(j,i)
2771 sigiso2(i,j)=sigiso2(j,i)
2772 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2773 nstate(i,j) = nstate(j,i)
2774 dtail(1,i,j) = dtail(1,j,i)
2775 dtail(2,i,j) = dtail(2,j,i)
2777 alphasur(k,i,j) = alphasur(k,j,i)
2778 wstate(k,i,j) = wstate(k,j,i)
2779 alphiso(k,i,j) = alphiso(k,j,i)
2782 dhead(2,1,i,j) = dhead(1,1,j,i)
2783 dhead(2,2,i,j) = dhead(1,2,j,i)
2784 dhead(1,1,i,j) = dhead(2,1,j,i)
2785 dhead(1,2,i,j) = dhead(2,2,j,i)
2787 epshead(i,j) = epshead(j,i)
2788 sig0head(i,j) = sig0head(j,i)
2791 wqdip(k,i,j) = wqdip(k,j,i)
2794 wquad(i,j) = wquad(j,i)
2795 epsintab(i,j) = epsintab(j,i)
2796 debaykap(i,j)=debaykap(j,i)
2797 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2804 !--------------------- GBV potential -----------------------------------
2806 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2807 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2808 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2809 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2811 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2812 write (iout,'(a/)') 'The epsilon array:'
2813 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2814 write (iout,'(/a)') 'One-body parameters:'
2815 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2816 's||/s_|_^2',' chip ',' alph '
2817 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2818 sigii(i),chip(i),alp(i),i=1,ntyp)
2821 write(iout,*)"Wrong ipot"
2827 !-----------------------------------------------------------------------
2828 ! Calculate the "working" parameters of SC interactions.
2830 !el from module energy - COMMON.INTERACT-------
2831 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2832 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2833 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2834 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2835 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2836 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2838 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2844 if (scelemode.eq.0) then
2853 sc_aa_tube_par(:)=0.0d0
2854 sc_bb_tube_par(:)=0.0d0
2856 !--------------------------------
2861 epslip(i,j)=epslip(j,i)
2864 if (scelemode.eq.0) then
2867 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2868 sigma(j,i)=sigma(i,j)
2869 rs0(i,j)=dwa16*sigma(i,j)
2874 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2875 'Working parameters of the SC interactions:',&
2876 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2881 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2883 ! print *,"SIGMA ZLA?",sigma(i,j)
2891 sigeps=dsign(1.0D0,epsij)
2893 aa_aq(i,j)=epsij*rrij*rrij
2894 ! print *,"ADASKO",epsij,rrij,expon
2895 bb_aq(i,j)=-sigeps*epsij*rrij
2896 aa_aq(j,i)=aa_aq(i,j)
2897 bb_aq(j,i)=bb_aq(i,j)
2898 epsijlip=epslip(i,j)
2899 sigeps=dsign(1.0D0,epsijlip)
2900 epsijlip=dabs(epsijlip)
2901 aa_lip(i,j)=epsijlip*rrij*rrij
2902 bb_lip(i,j)=-sigeps*epsijlip*rrij
2903 aa_lip(j,i)=aa_lip(i,j)
2904 bb_lip(j,i)=bb_lip(i,j)
2905 !C write(iout,*) aa_lip
2906 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2907 sigt1sq=sigma0(i)**2
2908 sigt2sq=sigma0(j)**2
2911 ratsig1=sigt2sq/sigt1sq
2912 ratsig2=1.0D0/ratsig1
2913 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2914 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2915 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2919 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2920 sigmaii(i,j)=rsum_max
2921 sigmaii(j,i)=rsum_max
2923 ! sigmaii(i,j)=r0(i,j)
2924 ! sigmaii(j,i)=r0(i,j)
2926 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2927 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2928 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2929 augm(i,j)=epsij*r_augm**(2*expon)
2930 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2937 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2938 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2939 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2944 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2945 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2946 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2947 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2948 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2949 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2950 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2951 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2952 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2953 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2954 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2955 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2956 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2957 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2958 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2959 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2968 read (isidep_nucl,*) ipot_nucl
2969 ! print *,"TU?!",ipot_nucl
2970 if (ipot_nucl.eq.1) then
2971 do i=1,ntyp_molec(2)
2972 do j=i,ntyp_molec(2)
2973 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2974 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2978 do i=1,ntyp_molec(2)
2979 do j=i,ntyp_molec(2)
2980 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2981 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2982 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2986 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2987 do i=1,ntyp_molec(2)
2988 do j=i,ntyp_molec(2)
2989 rrij=sigma_nucl(i,j)
2993 epsij=4*eps_nucl(i,j)
2994 sigeps=dsign(1.0D0,epsij)
2996 aa_nucl(i,j)=epsij*rrij*rrij
2997 bb_nucl(i,j)=-sigeps*epsij*rrij
2998 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2999 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
3000 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
3001 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
3002 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
3003 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
3006 aa_nucl(i,j)=aa_nucl(j,i)
3007 bb_nucl(i,j)=bb_nucl(j,i)
3008 ael3_nucl(i,j)=ael3_nucl(j,i)
3009 ael6_nucl(i,j)=ael6_nucl(j,i)
3010 ael63_nucl(i,j)=ael63_nucl(j,i)
3011 ael32_nucl(i,j)=ael32_nucl(j,i)
3012 elpp3_nucl(i,j)=elpp3_nucl(j,i)
3013 elpp6_nucl(i,j)=elpp6_nucl(j,i)
3014 elpp63_nucl(i,j)=elpp63_nucl(j,i)
3015 elpp32_nucl(i,j)=elpp32_nucl(j,i)
3016 eps_nucl(i,j)=eps_nucl(j,i)
3017 sigma_nucl(i,j)=sigma_nucl(j,i)
3018 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
3022 write(iout,*) "tube param"
3023 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
3024 ccavtubpep,dcavtubpep,tubetranenepep
3025 sigmapeptube=sigmapeptube**6
3026 sigeps=dsign(1.0D0,epspeptube)
3027 epspeptube=dabs(epspeptube)
3028 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
3029 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
3030 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
3032 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
3033 ccavtub(i),dcavtub(i),tubetranene(i)
3034 sigmasctube=sigmasctube**6
3035 sigeps=dsign(1.0D0,epssctube)
3036 epssctube=dabs(epssctube)
3037 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
3038 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
3039 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
3041 !-----------------READING SC BASE POTENTIALS-----------------------------
3042 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3043 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3044 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3045 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3046 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3047 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3048 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3049 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3050 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3051 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3052 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3053 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3054 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3055 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3056 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3057 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3059 write (iout,*) "ESCBASEPARM"
3060 do i=1,ntyp_molec(1)
3061 do j=1,ntyp_molec(2) ! without U then we will take T for U
3062 ! write (*,*) "Im in ", i, " ", j
3063 read(isidep_scbase,*) &
3064 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3065 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3066 ! write(*,*) "eps",eps_scbase(i,j)
3067 read(isidep_scbase,*) &
3068 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3069 chis_scbase(i,j,1),chis_scbase(i,j,2)
3070 read(isidep_scbase,*) &
3071 dhead_scbasei(i,j), &
3072 dhead_scbasej(i,j), &
3073 rborn_scbasei(i,j),rborn_scbasej(i,j)
3074 read(isidep_scbase,*) &
3075 (wdipdip_scbase(k,i,j),k=1,3), &
3076 (wqdip_scbase(k,i,j),k=1,2)
3077 read(isidep_scbase,*) &
3078 alphapol_scbase(i,j), &
3079 epsintab_scbase(i,j)
3080 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3081 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3082 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3083 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3084 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3085 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3086 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3087 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3089 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3090 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3091 write(*,*) "eps",eps_scbase(i,j)
3093 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3094 chis_scbase(i,j,1),chis_scbase(i,j,2)
3096 dhead_scbasei(i,j), &
3097 dhead_scbasej(i,j), &
3098 rborn_scbasei(i,j),rborn_scbasej(i,j)
3100 (wdipdip_scbase(k,i,j),k=1,3), &
3101 (wqdip_scbase(k,i,j),k=1,2)
3103 alphapol_scbase(i,j), &
3104 epsintab_scbase(i,j)
3109 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3110 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3111 write(*,*) "eps",eps_scbase(i,j)
3113 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3114 chis_scbase(i,j,1),chis_scbase(i,j,2)
3116 dhead_scbasei(i,j), &
3117 dhead_scbasej(i,j), &
3118 rborn_scbasei(i,j),rborn_scbasej(i,j)
3120 (wdipdip_scbase(k,i,j),k=1,3), &
3121 (wqdip_scbase(k,i,j),k=1,2)
3123 alphapol_scbase(i,j), &
3124 epsintab_scbase(i,j)
3127 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3128 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3130 do i=1,ntyp_molec(1)
3131 do j=1,ntyp_molec(2)-1
3132 epsij=eps_scbase(i,j)
3133 rrij=sigma_scbase(i,j)
3138 sigeps=dsign(1.0D0,epsij)
3140 aa_scbase(i,j)=epsij*rrij*rrij
3141 bb_scbase(i,j)=-sigeps*epsij*rrij
3144 !-----------------READING PEP BASE POTENTIALS-------------------
3145 allocate(eps_pepbase(ntyp_molec(2)))
3146 allocate(sigma_pepbase(ntyp_molec(2)))
3147 allocate(chi_pepbase(ntyp_molec(2),2))
3148 allocate(chipp_pepbase(ntyp_molec(2),2))
3149 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3150 allocate(sigmap1_pepbase(ntyp_molec(2)))
3151 allocate(sigmap2_pepbase(ntyp_molec(2)))
3152 allocate(chis_pepbase(ntyp_molec(2),2))
3153 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3155 write (iout,*) "EPEPBASEPARM"
3157 do j=1,ntyp_molec(2) ! without U then we will take T for U
3158 write (*,*) "Im in ", i, " ", j
3159 read(isidep_pepbase,*) &
3160 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3161 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3162 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3163 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3164 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3165 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3166 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3167 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3168 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3169 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3171 write(*,*) "eps",eps_pepbase(j)
3172 read(isidep_pepbase,*) &
3173 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3174 chis_pepbase(j,1),chis_pepbase(j,2)
3175 read(isidep_pepbase,*) &
3176 (wdipdip_pepbase(k,j),k=1,3)
3177 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3178 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3180 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3181 chis_pepbase(j,1),chis_pepbase(j,2)
3183 (wdipdip_pepbase(k,j),k=1,3)
3187 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3188 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3190 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3191 chis_pepbase(j,1),chis_pepbase(j,2)
3193 (wdipdip_pepbase(k,j),k=1,3)
3195 allocate(aa_pepbase(ntyp_molec(2)))
3196 allocate(bb_pepbase(ntyp_molec(2)))
3198 do j=1,ntyp_molec(2)-1
3199 epsij=eps_pepbase(j)
3200 rrij=sigma_pepbase(j)
3205 sigeps=dsign(1.0D0,epsij)
3207 aa_pepbase(j)=epsij*rrij*rrij
3208 bb_pepbase(j)=-sigeps*epsij*rrij
3210 !--------------READING SC PHOSPHATE-------------------------------------
3211 allocate(eps_scpho(ntyp_molec(1)))
3212 allocate(sigma_scpho(ntyp_molec(1)))
3213 allocate(chi_scpho(ntyp_molec(1),2))
3214 allocate(chipp_scpho(ntyp_molec(1),2))
3215 allocate(alphasur_scpho(4,ntyp_molec(1)))
3216 allocate(sigmap1_scpho(ntyp_molec(1)))
3217 allocate(sigmap2_scpho(ntyp_molec(1)))
3218 allocate(chis_scpho(ntyp_molec(1),2))
3219 allocate(wqq_scpho(ntyp_molec(1)))
3220 allocate(wqdip_scpho(2,ntyp_molec(1)))
3221 allocate(alphapol_scpho(ntyp_molec(1)))
3222 allocate(epsintab_scpho(ntyp_molec(1)))
3223 allocate(dhead_scphoi(ntyp_molec(1)))
3224 allocate(rborn_scphoi(ntyp_molec(1)))
3225 allocate(rborn_scphoj(ntyp_molec(1)))
3226 allocate(alphi_scpho(ntyp_molec(1)))
3230 do j=1,ntyp_molec(1) ! without U then we will take T for U
3231 write (*,*) "Im in scpho ", i, " ", j
3232 read(isidep_scpho,*) &
3233 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3234 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3235 write(*,*) "eps",eps_scpho(j)
3236 read(isidep_scpho,*) &
3237 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3238 chis_scpho(j,1),chis_scpho(j,2)
3239 read(isidep_scpho,*) &
3240 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3241 read(isidep_scpho,*) &
3242 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3244 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3245 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3246 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3247 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3248 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3249 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3250 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3251 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3255 allocate(aa_scpho(ntyp_molec(1)))
3256 allocate(bb_scpho(ntyp_molec(1)))
3258 do j=1,ntyp_molec(1)
3265 sigeps=dsign(1.0D0,epsij)
3267 aa_scpho(j)=epsij*rrij*rrij
3268 bb_scpho(j)=-sigeps*epsij*rrij
3272 read(isidep_peppho,*) &
3273 eps_peppho,sigma_peppho
3274 read(isidep_peppho,*) &
3275 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3276 read(isidep_peppho,*) &
3277 (wqdip_peppho(k),k=1,2)
3285 sigeps=dsign(1.0D0,epsij)
3287 aa_peppho=epsij*rrij*rrij
3288 bb_peppho=-sigeps*epsij*rrij
3291 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3296 ! Define the SC-p interaction constants (hard-coded; old style)
3299 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3301 ! aad(i,1)=0.3D0*4.0D0**12
3302 ! Following line for constants currently implemented
3303 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3304 aad(i,1)=1.5D0*4.0D0**12
3305 ! aad(i,1)=0.17D0*5.6D0**12
3307 ! "Soft" SC-p repulsion
3309 ! Following line for constants currently implemented
3310 ! aad(i,1)=0.3D0*4.0D0**6
3311 ! "Hard" SC-p repulsion
3312 bad(i,1)=3.0D0*4.0D0**6
3313 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3322 ! 8/9/01 Read the SC-p interaction constants from file
3325 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3328 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3329 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3330 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3331 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3335 write (iout,*) "Parameters of SC-p interactions:"
3337 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3338 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3343 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3345 do i=1,ntyp_molec(2)
3346 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3348 do i=1,ntyp_molec(2)
3349 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3350 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3352 r0pp=1.12246204830937298142*5.16158
3358 ! Define the constants of the disulfide bridge
3362 ! Old arbitrary potential - commented out.
3367 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3368 ! energy surface of diethyl disulfide.
3369 ! A. Liwo and U. Kozlowska, 11/24/03
3385 allocate(ichargelipid(ntyp_molec(4)))
3386 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3387 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3388 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3389 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3390 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3391 kjtokcal=0.2390057361
3393 !HERE THE MASS of MARTINI
3394 write(*,*) "before MARTINI PARAM"
3395 do i=1,ntyp_molec(4)
3401 !relative dielectric constant = 15 for implicit screening
3402 k_coulomb_lip=332.0d0/15.0d0
3403 !kbond = 1250 kJ/(mol*nm*2)
3404 kbondlip=1250.0d0*kjtokcal/100.0d0
3406 lip_angle_force=0.0d0
3407 if (DRY_MARTINI.gt.0) then
3408 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3409 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3410 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3411 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3412 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3413 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3414 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3415 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3417 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3418 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3419 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3420 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3421 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3422 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3423 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3424 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3426 lip_angle_angle=0.0d0
3427 lip_angle_angle(3,12,12)=120.0/krad
3428 lip_angle_angle(3,12,18)=180.0/krad
3429 lip_angle_angle(3,18,16)=180.0/krad
3430 lip_angle_angle(12,18,16)=180.0/krad
3431 lip_angle_angle(18,16,18)=120.0/krad
3432 lip_angle_angle(16,18,18)=180.0/krad
3433 lip_angle_angle(12,18,18)=180.0/krad
3434 lip_angle_angle(18,18,18)=180.0/krad
3435 read(ilipbond,*) temp1
3437 read(ilipbond,*) temp1, lip_bond(i,1), &
3438 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3439 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3440 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3441 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3444 lip_bond(i,j)=lip_bond(i,j)*10
3448 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3449 read(ilipnonbond,*) temp1
3451 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3452 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3453 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3454 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3455 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3457 ! write(*,*) i, lip_eps(i,18)
3459 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3462 read(ilipnonbond,*) temp1
3464 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3465 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3466 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3467 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3468 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3471 lip_sig(i,j)=lip_sig(i,j)*10.0
3474 write(*,*) "after MARTINI PARAM"
3478 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3479 allocate(alphapolcat2(ntyp,ntyp))
3480 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3481 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3482 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3483 allocate(epsintabcat(ntyp,ntyp))
3484 allocate(dtailcat(2,ntyp,ntyp))
3485 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3486 allocate(wqdipcat(2,ntyp,ntyp))
3487 allocate(wstatecat(4,ntyp,ntyp))
3488 allocate(dheadcat(2,2,ntyp,ntyp))
3489 allocate(nstatecat(ntyp,ntyp))
3490 allocate(debaykapcat(ntyp,ntyp))
3492 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3493 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3494 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3495 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3496 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3499 if (.not.allocated(ichargecat))&
3500 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
3501 write(*,*) "before ions",oldion
3504 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3505 if (oldion.eq.0) then
3506 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3507 allocate(icharge(1:ntyp1))
3508 read(iion,*) (icharge(i),i=1,ntyp)
3512 print *,ntyp_molec(5)
3513 do i=-ntyp_molec(5),ntyp_molec(5)
3514 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3515 print *,msc(i,5),restok(i,5)
3521 do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3523 ! do j=1,ntyp_molec(5)
3524 ! write (*,*) "Im in ALAB", i, " ", j
3526 epscat(i,j),sigmacat(i,j), &
3527 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3528 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3530 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3531 ! chiscat(i,j),chiscat(j,i), &
3532 chis1cat(i,j),chis2cat(i,j), &
3534 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3535 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3536 dtailcat(1,i,j),dtailcat(2,i,j), &
3537 epsheadcat(i,j),sig0headcat(i,j), &
3539 ! rborncat(i,j),rborncat(j,i),&
3540 rborn1cat(i,j),rborn2cat(i,j),&
3541 (wqdipcat(k,i,j),k=1,2), &
3542 alphapolcat(i,j),alphapolcat2(j,i), &
3543 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3545 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3546 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3548 ! write (iout,*) 'i= ', i, ' j= ', j
3549 ! write (iout,*) 'epsi0= ', epscat(1,j)
3550 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3551 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3552 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3553 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3554 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3555 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3556 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3557 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3558 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3559 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3560 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3561 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3562 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3563 ! write (iout,*) 'a1= ', rborncat(j,1)
3564 ! write (iout,*) 'a2= ', rborncat(1,j)
3565 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3566 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3567 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3568 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3572 ! If ((i.eq.1).and.(j.eq.27)) then
3573 ! write (iout,*) 'SEP'
3574 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3575 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3580 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3582 do j=1,ntyp_molec(5)
3586 sigeps=dsign(1.0D0,epsij)
3588 aa_aq_cat(i,j)=epsij*rrij*rrij
3589 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3594 do j=1,ntyp_molec(5)-1
3596 write (iout,*) 'i= ', i, ' j= ', j
3597 write (iout,*) 'epsi0= ', epscat(i,j)
3598 write (iout,*) 'sigma0= ', sigmacat(i,j)
3599 write (iout,*) 'chi1= ', chi1cat(i,j)
3600 write (iout,*) 'chi1= ', chi2cat(i,j)
3601 write (iout,*) 'chip1= ', chipp1cat(i,j)
3602 write (iout,*) 'chip2= ', chipp2cat(i,j)
3603 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3604 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3605 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3606 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3607 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3608 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3609 write (iout,*) 'chis1= ', chis1cat(i,j)
3610 write (iout,*) 'chis1= ', chis2cat(i,j)
3611 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3612 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3613 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3614 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3615 write (iout,*) 'a1= ', rborn1cat(i,j)
3616 write (iout,*) 'a2= ', rborn2cat(i,j)
3617 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3618 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3619 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3620 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3621 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3622 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3625 If ((i.eq.1).and.(j.eq.27)) then
3626 write (iout,*) 'SEP'
3627 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3628 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3635 ! read number of Zn2+
3636 ! here two denotes the Zn2+ and Cu2+
3637 write(iout,*) "before TRANPARM"
3638 allocate(aomicattr(0:3,2))
3639 allocate(athetacattran(0:6,5,2))
3640 allocate(agamacattran(3,5,2))
3641 allocate(acatshiftdsc(5,2))
3642 allocate(bcatshiftdsc(5,2))
3643 allocate(demorsecat(5,2))
3644 allocate(alphamorsecat(5,2))
3645 allocate(x0catleft(5,2))
3646 allocate(x0catright(5,2))
3647 allocate(x0cattrans(5,2))
3648 allocate(ntrantyp(2))
3649 do i=1,1 ! currently only Zn2+
3651 read(iiontran,*) ntrantyp(i)
3654 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3655 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3656 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3657 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3658 read(iiontran,*) (aomicattr(j,i),j=0,3)
3660 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3661 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3662 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3667 write (iout,'(/a)') "Disulfide bridge parameters:"
3668 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3669 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3670 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3671 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3674 if (shield_mode.gt.0) then
3675 pi=4.0D0*datan(1.0D0)
3676 !C VSolvSphere the volume of solving sphere
3678 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3679 !C there will be no distinction between proline peptide group and normal peptide
3680 !C group in case of shielding parameters
3681 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3682 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3683 write (iout,*) VSolvSphere,VSolvSphere_div
3684 !C long axis of side chain
3686 long_r_sidechain(i)=vbldsc0(1,i)
3687 ! if (scelemode.eq.0) then
3688 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3689 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3691 ! short_r_sidechain(i)=sigma(i,i)
3693 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3700 111 write (iout,*) "Error reading bending energy parameters."
3702 112 write (iout,*) "Error reading rotamer energy parameters."
3704 113 write (iout,*) "Error reading torsional energy parameters."
3706 114 write (iout,*) "Error reading double torsional energy parameters."
3708 115 write (iout,*) &
3709 "Error reading cumulant (multibody energy) parameters."
3711 116 write (iout,*) "Error reading electrostatic energy parameters."
3713 117 write (iout,*) "Error reading side chain interaction parameters."
3715 118 write (iout,*) "Error reading SCp interaction parameters."
3717 119 write (iout,*) "Error reading SCCOR parameters"
3719 121 write (iout,*) "Error in Czybyshev parameters"
3721 123 write(iout,*) "Error in transition metal parameters"
3724 call MPI_Finalize(Ierror)
3728 end subroutine parmread
3730 !-----------------------------------------------------------------------------
3732 !-----------------------------------------------------------------------------
3733 subroutine printmat(ldim,m,n,iout,key,a)
3736 character(len=3),dimension(n) :: key
3737 real(kind=8),dimension(ldim,n) :: a
3739 integer :: i,j,k,m,iout,nlim
3743 write (iout,1000) (key(k),k=i,nlim)
3745 1000 format (/5x,8(6x,a3))
3746 1020 format (/80(1h-)/)
3748 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3751 1010 format (a3,2x,8(f9.4))
3753 end subroutine printmat
3754 !-----------------------------------------------------------------------------
3756 !-----------------------------------------------------------------------------
3758 ! Read the PDB file and convert the peptide geometry into virtual-chain
3761 use energy_data, only: itype,istype
3765 ! use control, only: rescode,sugarcode
3766 ! implicit real*8 (a-h,o-z)
3767 ! include 'DIMENSIONS'
3768 ! include 'COMMON.LOCAL'
3769 ! include 'COMMON.VAR'
3770 ! include 'COMMON.CHAIN'
3771 ! include 'COMMON.INTERACT'
3772 ! include 'COMMON.IOUNITS'
3773 ! include 'COMMON.GEO'
3774 ! include 'COMMON.NAMES'
3775 ! include 'COMMON.CONTROL'
3776 ! include 'COMMON.DISTFIT'
3777 ! include 'COMMON.SETUP'
3778 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3780 logical :: lprn=.true.,fail
3781 real(kind=8),dimension(3) :: e1,e2,e3
3782 real(kind=8) :: dcj,efree_temp
3783 character(len=3) :: seq,res,res2
3784 character(len=5) :: atom
3785 character(len=80) :: card
3786 real(kind=8),dimension(3,20) :: sccor
3787 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3788 integer :: isugar,molecprev,firstion
3789 character*1 :: sugar
3791 real(kind=8),dimension(3) :: ccc
3793 integer,dimension(2,maxres/3) :: hfrag_alloc
3794 integer,dimension(4,maxres/3) :: bfrag_alloc
3795 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3796 real(kind=8),dimension(:,:), allocatable :: c_temporary
3797 integer,dimension(:,:) , allocatable :: itype_temporary
3798 integer,dimension(:),allocatable :: istype_temp
3805 ! write (2,*) "UNRES_PDB",unres_pdb
3825 !-----------------------------
3826 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3827 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3828 if(.not. allocated(istype)) allocate(istype(maxres))
3830 read (ipdbin,'(a80)',end=10) card
3831 write (iout,'(a)') card
3832 if (card(:5).eq.'HELIX') then
3835 read(card(22:25),*) hfrag(1,nhfrag)
3836 read(card(34:37),*) hfrag(2,nhfrag)
3838 if (card(:5).eq.'SHEET') then
3841 read(card(24:26),*) bfrag(1,nbfrag)
3842 read(card(35:37),*) bfrag(2,nbfrag)
3843 !rc----------------------------------------
3844 !rc to be corrected !!!
3845 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3846 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3847 !rc----------------------------------------
3849 if (card(:3).eq.'END') then
3851 else if (card(:3).eq.'TER') then
3856 itype(ires_old,molecule)=ntyp1_molec(molecule)
3857 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3858 nres_molec(molecule)=nres_molec(molecule)+2
3859 ! if (molecule.eq.4) ires=ires+2
3861 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3864 dc(j,ires)=sccor(j,iii)
3867 call sccenter(ires,iii,sccor)
3871 else if (card(:3).eq.'BRA') then
3875 itype(ires,molecule)=ntyp1_molec(molecule)-1
3876 nres_molec(molecule)=nres_molec(molecule)+1
3880 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3881 ! Fish out the ATOM cards.
3882 ! write(iout,*) 'card',card(1:20)
3883 ! print *,"ATU ",card(1:6), CARD(3:6)
3885 if (index(card(1:4),'ATOM').gt.0) then
3886 read (card(12:16),*) atom
3887 ! write (iout,*) "! ",atom," !",ires
3888 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3889 if (card(14:16).eq.'LIP') then
3894 nres_molec(molecule)=nres_molec(molecule)+1
3895 itype(ires,molecule)=ntyp1_molec(molecule)
3904 nres_molec(molecule)=nres_molec(molecule)+1
3905 read (card(18:20),'(a3)') res
3907 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3908 if (UNRES_PDB) then!
3909 itype(ires,molecule)=rescode(ires,res,0,molecule)
3911 itype(ires,molecule)=rescode_lip(res,ires)
3914 read (card(23:26),*) ires
3915 read (card(18:20),'(a3)') res
3916 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3917 ! & " ires_old",ires_old
3918 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3919 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3920 if (ires-ishift+ishift1.ne.ires_old) then
3921 ! Calculate the CM of the preceding residue.
3922 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3924 ! write (iout,*) "Calculating sidechain center iii",iii
3927 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3930 call sccenter(ires_old,iii,sccor)
3934 ! Start new residue.
3935 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3938 else if (ibeg.eq.1) then
3939 write (iout,*) "BEG ires",ires
3941 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3944 nres_molec(molecule)=nres_molec(molecule)+1
3946 ires=ires-ishift+ishift1
3948 ! write (iout,*) "ishift",ishift," ires",ires,&
3949 ! " ires_old",ires_old
3951 else if (ibeg.eq.2) then
3953 ishift=-ires_old+ires-1 !!!!!
3954 ishift1=ishift1-1 !!!!!
3955 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3956 ires=ires-ishift+ishift1
3957 ! print *,ires,ishift,ishift1
3961 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3962 ires=ires-ishift+ishift1
3965 ! print *,'atom',ires,atom
3966 if (res.eq.'ACE' .or. res.eq.'NHE') then
3969 if (atom.eq.'CA '.or.atom.eq.'N ') then
3971 itype(ires,molecule)=rescode(ires,res,0,molecule)
3973 ! nres_molec(molecule)=nres_molec(molecule)+1
3977 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3978 ! nres_molec(molecule)=nres_molec(molecule)+1
3979 read (card(19:19),'(a1)') sugar
3980 isugar=sugarcode(sugar,ires)
3981 ! if (ibeg.eq.1) then
3985 ! print *,"ires=",ires,istype(ires)
3991 ires=ires-ishift+ishift1
3993 ! write (iout,*) "ires_old",ires_old," ires",ires
3994 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3997 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3998 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3999 res.eq.'NHE'.and.atom(:2).eq.'HN') then
4000 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4001 ! print *,ires,ishift,ishift1
4002 ! write (iout,*) "backbone ",atom
4004 write (iout,'(2i3,2x,a,3f8.3)') &
4005 ires,itype(ires,1),res,(c(j,ires),j=1,3)
4008 nres_molec(molecule)=nres_molec(molecule)+1
4010 sccor(j,iii)=c(j,ires)
4012 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
4013 atom.eq."C2'" .or. atom.eq."C3'" &
4014 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
4015 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
4016 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
4017 ! print *,ires,ishift,ishift1
4021 ! sccor(j,iii)=c(j,ires)
4024 c(j,ires)=c(j,ires)+ccc(j)/5.0
4026 print *,counter,molecule
4027 if (counter.eq.5) then
4029 nres_molec(molecule)=nres_molec(molecule)+1
4032 ! sccor(j,iii)=c(j,ires)
4036 ! print *, "ATOM",atom(1:3)
4037 else if (atom.eq."C5'") then
4038 read (card(19:19),'(a1)') sugar
4039 isugar=sugarcode(sugar,ires)
4044 ! print *,ires,istype(ires)
4048 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4049 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4050 nres_molec(molecule)=nres_molec(molecule)+1
4051 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4055 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4057 else if ((atom.eq."C1'").and.unres_pdb) then
4059 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4060 ! write (*,*) card(23:27),ires,itype(ires,1)
4061 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4062 atom.ne.'N' .and. atom.ne.'C' .and. &
4063 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4064 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4065 .and. atom.ne.'P '.and. &
4066 atom(1:1).ne.'H' .and. &
4067 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4069 ! write (iout,*) "sidechain ",atom
4070 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4071 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4072 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4074 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4078 ! print *,"IONS",ions,card(1:6)
4079 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4080 if (firstion.eq.0) then
4084 dc(j,ires)=sccor(j,iii)
4087 call sccenter(ires,iii,sccor)
4090 read (card(12:16),*) atom
4091 print *,"HETATOM", atom(1:2)
4092 read (card(18:20),'(a3)') res
4093 if (atom(3:3).eq.'H') cycle
4094 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4095 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4096 .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
4097 (atom(1:2).eq.'O ')) then
4099 print *,"I have water"
4100 if (molecule.ne.5) molecprev=molecule
4102 nres_molec(molecule)=nres_molec(molecule)+1
4103 print *,"HERE",nres_molec(molecule)
4104 if (res.ne.'WAT') res=res(2:3)//' '
4105 itype(ires,molecule)=rescode(ires,res,0,molecule)
4106 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4110 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4111 if (ires.eq.0) return
4112 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4115 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4117 nres_molec(molecule)=nres_molec(molecule)-2
4118 print *,'I have',nres, nres_molec(:)
4120 do k=1,4 ! ions are without dummy
4121 if (nres_molec(k).eq.0) cycle
4123 ! write (iout,*) i,itype(i,1)
4124 ! if (itype(i,1).eq.ntyp1) then
4125 ! write (iout,*) "dummy",i,itype(i,1)
4127 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4128 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4132 if (itype(i,k).eq.ntyp1_molec(k)) then
4133 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4134 if (itype(i+2,k).eq.0) then
4135 ! print *,"masz sieczke"
4137 if (itype(i+2,j).ne.0) then
4139 itype(i+1,j)=ntyp1_molec(j)
4140 nres_molec(k)=nres_molec(k)-1
4141 nres_molec(j)=nres_molec(j)+1
4147 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4148 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4149 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4150 ! if (unres_pdb) then
4151 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4152 ! print *,i,'tu dochodze'
4153 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4161 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4165 dcj=(c(j,i-2)-c(j,i-3))/2.0
4166 if (dcj.eq.0) dcj=1.23591524223
4171 else !itype(i+1,1).eq.ntyp1
4172 ! if (unres_pdb) then
4173 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4174 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4181 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4182 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4186 dcj=(c(j,i+3)-c(j,i+2))/2.0
4187 if (dcj.eq.0) dcj=1.23591524223
4192 endif !itype(i+1,1).eq.ntyp1
4193 endif !itype.eq.ntyp1
4197 ! Calculate the CM of the last side chain.
4201 dc(j,ires)=sccor(j,iii)
4204 call sccenter(ires,iii,sccor)
4210 ! print *,"molecule",molecule
4211 if ((itype(nres,1).ne.10)) then
4214 if (molecule.eq.5) molecule=molecprev
4215 itype(nres,molecule)=ntyp1_molec(molecule)
4216 nres_molec(molecule)=nres_molec(molecule)+1
4217 ! if (unres_pdb) then
4218 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4219 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4226 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4230 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4231 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4232 c(j,2*nres)=c(j,nres)
4236 ! print *,'I have',nres, nres_molec(:)
4238 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4240 if (nres.ne.nres0) then
4241 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4243 stop "Error nres value in WHAM input"
4246 !---------------------------------
4247 !el reallocate tables
4250 ! hfrag_alloc(j,i)=hfrag(j,i)
4253 ! bfrag_alloc(j,i)=bfrag(j,i)
4259 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4260 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4261 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4262 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4266 ! hfrag(j,i)=hfrag_alloc(j,i)
4271 ! bfrag(j,i)=bfrag_alloc(j,i)
4274 !el end reallocate tables
4275 !---------------------------------
4283 c(j,2*nres)=c(j,nres)
4286 if (itype(1,1).eq.ntyp1) then
4290 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4291 call refsys(2,3,4,e1,e2,e3,fail)
4298 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4299 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4303 dcj=(c(j,4)-c(j,3))/2.0
4309 ! First lets assign correct dummy to correct type of chain
4311 if (itype(1,1).eq.ntyp1) then
4312 if (itype(2,1).eq.0) then
4314 if (itype(2,j).ne.0) then
4316 itype(1,j)=ntyp1_molec(j)
4317 nres_molec(1)=nres_molec(1)-1
4318 nres_molec(j)=nres_molec(j)+1
4325 print *,'I have',nres, nres_molec(:)
4327 ! Copy the coordinates to reference coordinates
4333 ! Calculate internal coordinates.
4335 write (iout,'(/a)') &
4336 "Cartesian coordinates of the reference structure"
4337 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4338 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4340 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4341 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4342 (c(j,ires+nres),j=1,3)
4345 ! znamy już nres wiec mozna alokowac tablice
4346 ! Calculate internal coordinates.
4347 if(me.eq.king.or..not.out1file)then
4348 write (iout,'(a)') &
4349 "Backbone and SC coordinates as read from the PDB"
4351 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4352 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4353 (c(j,nres+ires),j=1,3)
4356 ! NOW LETS ROCK! SORTING
4357 allocate(c_temporary(3,2*nres))
4358 allocate(itype_temporary(nres,5))
4359 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4360 if (.not.allocated(istype)) write(iout,*) &
4361 "SOMETHING WRONG WITH ISTYTPE"
4362 allocate(istype_temp(nres))
4363 itype_temporary(:,:)=0
4367 if (itype(i,k).ne.0) then
4369 c_temporary(j,seqalingbegin)=c(j,i)
4370 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4372 itype_temporary(seqalingbegin,k)=itype(i,k)
4373 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4374 istype_temp(seqalingbegin)=istype(i)
4375 molnum(seqalingbegin)=k
4376 seqalingbegin=seqalingbegin+1
4382 c(j,i)=c_temporary(j,i)
4384 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4390 itype(i,k)=itype_temporary(i,k)
4391 istype(i)=istype_temp(i)
4394 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4395 ! I have only ions now dummy atoms in the system
4397 itype(1,5)=itype(2,5)
4403 itype(i,5)=itype(i+1,5)
4410 nres_molec(1)=nres_molec(1)-1
4412 ! if (itype(1,1).eq.ntyp1) then
4415 ! if (unres_pdb) then
4416 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4417 ! call refsys(2,3,4,e1,e2,e3,fail)
4424 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4428 ! dcj=(c(j,4)-c(j,3))/2.0
4430 ! c(j,nres+1)=c(j,1)
4436 write (iout,'(/a)') &
4437 "Cartesian coordinates of the reference structure after sorting"
4438 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4439 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4441 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4442 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4443 (c(j,ires+nres),j=1,3)
4447 print *,seqalingbegin,nres
4448 if(.not.allocated(vbld)) then
4449 allocate(vbld(2*nres))
4454 if(.not.allocated(vbld_inv)) then
4455 allocate(vbld_inv(2*nres))
4461 if(.not.allocated(theta)) then
4462 allocate(theta(nres+2))
4466 if(.not.allocated(phi)) allocate(phi(nres+2))
4467 if(.not.allocated(alph)) allocate(alph(nres+2))
4468 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4469 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4470 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4471 if(.not.allocated(costtab)) allocate(costtab(nres))
4472 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4473 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4474 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4475 if(.not.allocated(xxref)) allocate(xxref(nres))
4476 if(.not.allocated(yyref)) allocate(yyref(nres))
4477 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4478 if(.not.allocated(dc_norm)) then
4479 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4480 allocate(dc_norm(3,0:2*nres+2))
4483 write(iout,*) "before int_from_cart"
4484 call int_from_cart(.true.,.false.)
4485 call sc_loc_geom(.false.)
4486 write(iout,*) "after int_from_cart"
4490 thetaref(i)=theta(i)
4493 write(iout,*) "after thetaref"
4501 dc(j,i)=c(j,i+1)-c(j,i)
4502 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4507 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4508 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4510 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4514 ! Copy the coordinates to reference coordinates
4515 ! Splits to single chain if occurs
4519 ! cref(j,i,cou)=c(j,i)
4525 ! chomo(j,i,k)=c(j,i)
4528 ! write(iout,*) "after chomo"
4530 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4531 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4532 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4533 !-----------------------------
4537 write (iout,*) "symetr", symetr
4540 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4542 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4545 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4550 cref(j,i,cou)=c(j,i)
4551 cref(j,i+nres,cou)=c(j,i+nres)
4553 chain_rep(j,lll,kkk)=c(j,i)
4554 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4558 write (iout,*) chain_length
4559 if (chain_length.eq.0) chain_length=nres
4561 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4562 chain_rep(j,chain_length+nres,symetr) &
4563 =chain_rep(j,chain_length+nres,1)
4566 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4568 ! do kkk=1,chain_length
4569 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4573 ! makes copy of chains
4574 write (iout,*) "symetr", symetr
4579 if (symetr.gt.1) then
4586 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4592 write (iout,*) i,icha
4593 do lll=1,chain_length
4595 if (cou.le.nres) then
4597 kupa=mod(lll,chain_length)
4598 iprzes=(kkk-1)*chain_length+lll
4599 if (kupa.eq.0) kupa=chain_length
4600 write (iout,*) "kupa", kupa
4601 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4602 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4609 !-koniec robienia kopii
4612 write (iout,*) "nowa struktura", nperm
4614 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4616 cref(3,i,kkk),cref(1,nres+i,kkk),&
4617 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4619 100 format (//' alpha-carbon coordinates ',&
4620 ' centroid coordinates'/ &
4621 ' ', 6X,'X',11X,'Y',11X,'Z', &
4622 10X,'X',11X,'Y',11X,'Z')
4623 110 format (a,'(',i5,')',6f12.5)
4629 bfrag(i,j)=bfrag(i,j)-ishift
4635 hfrag(i,j)=hfrag(i,j)-ishift
4641 end subroutine readpdb
4642 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4643 !-----------------------------------------------------------------------------
4645 !-----------------------------------------------------------------------------
4646 subroutine read_control
4660 use random, only: random_init
4661 ! implicit real*8 (a-h,o-z)
4662 ! include 'DIMENSIONS'
4664 use prng, only:prng_restart
4666 logical :: OKRandom!, prng_restart
4669 ! include 'COMMON.IOUNITS'
4670 ! include 'COMMON.TIME1'
4671 ! include 'COMMON.THREAD'
4672 ! include 'COMMON.SBRIDGE'
4673 ! include 'COMMON.CONTROL'
4674 ! include 'COMMON.MCM'
4675 ! include 'COMMON.MAP'
4676 ! include 'COMMON.HEADER'
4677 ! include 'COMMON.CSA'
4678 ! include 'COMMON.CHAIN'
4679 ! include 'COMMON.MUCA'
4680 ! include 'COMMON.MD'
4681 ! include 'COMMON.FFIELD'
4682 ! include 'COMMON.INTERACT'
4683 ! include 'COMMON.SETUP'
4684 !el integer :: KDIAG,ICORFL,IXDR
4685 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4686 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4687 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4688 ! character(len=80) :: ucase
4689 character(len=640) :: controlcard
4691 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4693 usampl=.false. ! the default value of usample should be 0
4694 ! write(iout,*) "KURWA2",usampl
4698 read (INP,'(a)') titel
4699 call card_concat(controlcard,.true.)
4700 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4701 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4702 call reada(controlcard,'SEED',seed,0.0D0)
4703 call random_init(seed)
4704 ! Set up the time limit (caution! The time must be input in minutes!)
4705 read_cart=index(controlcard,'READ_CART').gt.0
4706 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4707 call readi(controlcard,'SYM',symetr,1)
4708 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4709 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4710 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4711 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4712 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4713 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4714 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4715 call reada(controlcard,'DRMS',drms,0.1D0)
4716 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4717 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4718 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4719 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4720 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4721 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4722 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4723 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4724 write (iout,'(a,f10.1)')'DRMS = ',drms
4725 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4726 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4728 call readi(controlcard,'NZ_START',nz_start,0)
4729 call readi(controlcard,'NZ_END',nz_end,0)
4730 call readi(controlcard,'IZ_SC',iz_sc,0)
4731 timlim=60.0D0*timlim
4732 safety = 60.0d0*safety
4735 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4736 !C SHIELD keyword sets if the shielding effect of side-chains is used
4737 !C 0 denots no shielding is used all peptide are equally despite the
4738 !C solvent accesible area
4739 !C 1 the newly introduced function
4740 !C 2 reseved for further possible developement
4741 call readi(controlcard,'SHIELD',shield_mode,0)
4742 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4743 write(iout,*) "shield_mode",shield_mode
4744 call readi(controlcard,'TORMODE',tor_mode,0)
4745 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4746 write(iout,*) "torsional and valence angle mode",tor_mode
4748 !C Varibles set size of box
4749 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4750 protein=index(controlcard,"PROTEIN").gt.0
4751 ions=index(controlcard,"IONS").gt.0
4752 fodson=index(controlcard,"FODSON").gt.0
4753 call readi(controlcard,'OLDION',oldion,1)
4754 nucleic=index(controlcard,"NUCLEIC").gt.0
4755 write (iout,*) "with_theta_constr ",with_theta_constr
4756 AFMlog=(index(controlcard,'AFM'))
4757 selfguide=(index(controlcard,'SELFGUIDE'))
4758 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4759 call readi(controlcard,'GENCONSTR',genconstr,0)
4760 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4761 call reada(controlcard,'BOXY',boxysize,100.0d0)
4762 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4763 call readi(controlcard,'TUBEMOD',tubemode,0)
4764 print *,"SCELE",scelemode
4765 call readi(controlcard,"SCELEMODE",scelemode,0)
4766 print *,"SCELE",scelemode
4768 ! elemode = 0 is orignal UNRES electrostatics
4769 ! elemode = 1 is "Momo" potentials in progress
4770 ! elemode = 2 is in development EVALD
4773 write (iout,*) TUBEmode,"TUBEMODE"
4774 if (TUBEmode.gt.0) then
4775 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4776 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4777 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4778 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4779 call reada(controlcard,"VNANO",velnanoconst,0.0d0)
4780 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4781 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4782 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4783 buftubebot=bordtubebot+tubebufthick
4784 buftubetop=bordtubetop-tubebufthick
4787 ! CUTOFFF ON ELECTROSTATICS
4788 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4789 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4790 write(iout,*) "R_CUT_ELE=",r_cut_ele
4791 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4792 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4793 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4795 ! Lipidec parameters
4796 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4797 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4798 if (lipthick.gt.0.0d0) then
4799 bordliptop=(boxzsize+lipthick)/2.0
4800 bordlipbot=bordliptop-lipthick
4801 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4802 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4803 buflipbot=bordlipbot+lipbufthick
4804 bufliptop=bordliptop-lipbufthick
4805 if ((lipbufthick*2.0d0).gt.lipthick) &
4806 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4807 endif !lipthick.gt.0
4808 write(iout,*) "bordliptop=",bordliptop
4809 write(iout,*) "bordlipbot=",bordlipbot
4810 write(iout,*) "bufliptop=",bufliptop
4811 write(iout,*) "buflipbot=",buflipbot
4812 write (iout,*) "SHIELD MODE",shield_mode
4814 !C-------------------------
4815 minim=(index(controlcard,'MINIMIZE').gt.0)
4816 dccart=(index(controlcard,'CART').gt.0)
4817 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4818 overlapsc=.not.overlapsc
4819 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4820 searchsc=.not.searchsc
4821 sideadd=(index(controlcard,'SIDEADD').gt.0)
4822 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4823 outpdb=(index(controlcard,'PDBOUT').gt.0)
4824 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4825 pdbref=(index(controlcard,'PDBREF').gt.0)
4826 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4827 indpdb=index(controlcard,'PDBSTART')
4828 extconf=(index(controlcard,'EXTCONF').gt.0)
4829 call readi(controlcard,'IPRINT',iprint,0)
4830 call readi(controlcard,'MAXGEN',maxgen,10000)
4831 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4832 call readi(controlcard,"KDIAG",kdiag,0)
4833 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4834 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4835 write (iout,*) "RESCALE_MODE",rescale_mode
4836 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4837 if (index(controlcard,'REGULAR').gt.0.0D0) then
4838 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4842 if (index(controlcard,'CHECKGRAD').gt.0) then
4844 if (index(controlcard,'CART').gt.0) then
4846 elseif (index(controlcard,'CARINT').gt.0) then
4851 elseif (index(controlcard,'THREAD').gt.0) then
4853 call readi(controlcard,'THREAD',nthread,0)
4854 if (nthread.gt.0) then
4855 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4858 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4859 stop 'Error termination in Read_Control.'
4861 else if (index(controlcard,'MCMA').gt.0) then
4863 else if (index(controlcard,'MCEE').gt.0) then
4865 else if (index(controlcard,'MULTCONF').gt.0) then
4867 else if (index(controlcard,'MAP').gt.0) then
4869 call readi(controlcard,'MAP',nmap,0)
4870 else if (index(controlcard,'CSA').gt.0) then
4872 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4874 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4877 !fcm else if (index(controlcard,'MCMF').gt.0) then
4879 else if (index(controlcard,'SOFTREG').gt.0) then
4881 else if (index(controlcard,'CHECK_BOND').gt.0) then
4883 else if (index(controlcard,'TEST').gt.0) then
4885 else if (index(controlcard,'MD').gt.0) then
4887 else if (index(controlcard,'RE ').gt.0) then
4891 lmuca=index(controlcard,'MUCA').gt.0
4892 call readi(controlcard,'MUCADYN',mucadyn,0)
4893 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4894 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4896 write (iout,*) 'MUCADYN=',mucadyn
4897 write (iout,*) 'MUCASMOOTH=',muca_smooth
4900 iscode=index(controlcard,'ONE_LETTER')
4901 indphi=index(controlcard,'PHI')
4902 indback=index(controlcard,'BACK')
4903 iranconf=index(controlcard,'RAND_CONF')
4904 i2ndstr=index(controlcard,'USE_SEC_PRED')
4905 gradout=index(controlcard,'GRADOUT').gt.0
4906 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4907 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4908 if (me.eq.king .or. .not.out1file ) &
4909 write (iout,*) "DISTCHAINMAX",distchainmax
4911 if(me.eq.king.or..not.out1file) &
4912 write (iout,'(2a)') diagmeth(kdiag),&
4913 ' routine used to diagonalize matrices.'
4914 if (shield_mode.gt.0) then
4915 pi=4.0D0*datan(1.0D0)
4916 !C VSolvSphere the volume of solving sphere
4918 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4919 !C there will be no distinction between proline peptide group and normal peptide
4920 !C group in case of shielding parameters
4921 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4922 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4923 write (iout,*) VSolvSphere,VSolvSphere_div
4924 !C long axis of side chain
4926 ! long_r_sidechain(i)=vbldsc0(1,i)
4927 ! short_r_sidechain(i)=sigma0(i)
4928 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4935 end subroutine read_control
4936 !-----------------------------------------------------------------------------
4937 subroutine read_REMDpar
4939 ! Read REMD settings
4946 use control_data, only:out1file
4947 ! implicit real*8 (a-h,o-z)
4948 ! include 'DIMENSIONS'
4949 ! include 'COMMON.IOUNITS'
4950 ! include 'COMMON.TIME1'
4951 ! include 'COMMON.MD'
4954 !el include 'COMMON.LANGEVIN'
4956 !el include 'COMMON.LANGEVIN.lang0'
4958 ! include 'COMMON.INTERACT'
4959 ! include 'COMMON.NAMES'
4960 ! include 'COMMON.GEO'
4961 ! include 'COMMON.REMD'
4962 ! include 'COMMON.CONTROL'
4963 ! include 'COMMON.SETUP'
4964 ! character(len=80) :: ucase
4965 character(len=320) :: controlcard
4966 character(len=3200) :: controlcard1
4967 integer :: iremd_m_total
4970 ! real(kind=8) :: var,ene
4972 if(me.eq.king.or..not.out1file) &
4973 write (iout,*) "REMD setup"
4975 call card_concat(controlcard,.true.)
4976 call readi(controlcard,"NREP",nrep,3)
4977 call readi(controlcard,"NSTEX",nstex,1000)
4978 call reada(controlcard,"RETMIN",retmin,10.0d0)
4979 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4980 mremdsync=(index(controlcard,'SYNC').gt.0)
4981 call readi(controlcard,"NSYN",i_sync_step,100)
4982 restart1file=(index(controlcard,'REST1FILE').gt.0)
4983 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4984 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4985 if(max_cache_traj_use.gt.max_cache_traj) &
4986 max_cache_traj_use=max_cache_traj
4987 if(me.eq.king.or..not.out1file) then
4988 !d if (traj1file) then
4989 !rc caching is in testing - NTWX is not ignored
4990 !d write (iout,*) "NTWX value is ignored"
4991 !d write (iout,*) " trajectory is stored to one file by master"
4992 !d write (iout,*) " before exchange at NSTEX intervals"
4994 write (iout,*) "NREP= ",nrep
4995 write (iout,*) "NSTEX= ",nstex
4996 write (iout,*) "SYNC= ",mremdsync
4997 write (iout,*) "NSYN= ",i_sync_step
4998 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
5001 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
5002 if (index(controlcard,'TLIST').gt.0) then
5004 call card_concat(controlcard1,.true.)
5005 read(controlcard1,*) (remd_t(i),i=1,nrep)
5006 if(me.eq.king.or..not.out1file) &
5007 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
5010 if (index(controlcard,'MLIST').gt.0) then
5012 call card_concat(controlcard1,.true.)
5013 read(controlcard1,*) (remd_m(i),i=1,nrep)
5014 if(me.eq.king.or..not.out1file) then
5015 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
5018 iremd_m_total=iremd_m_total+remd_m(i)
5020 write (iout,*) 'Total number of replicas ',iremd_m_total
5023 if(me.eq.king.or..not.out1file) &
5024 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
5026 end subroutine read_REMDpar
5027 !-----------------------------------------------------------------------------
5028 subroutine read_MDpar
5032 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
5034 use geometry_data, only: pi
5036 ! implicit real*8 (a-h,o-z)
5037 ! include 'DIMENSIONS'
5038 ! include 'COMMON.IOUNITS'
5039 ! include 'COMMON.TIME1'
5040 ! include 'COMMON.MD'
5043 !el include 'COMMON.LANGEVIN'
5045 !el include 'COMMON.LANGEVIN.lang0'
5047 ! include 'COMMON.INTERACT'
5048 ! include 'COMMON.NAMES'
5049 ! include 'COMMON.GEO'
5050 ! include 'COMMON.SETUP'
5051 ! include 'COMMON.CONTROL'
5052 ! include 'COMMON.SPLITELE'
5053 ! character(len=80) :: ucase
5054 character(len=320) :: controlcard
5059 call card_concat(controlcard,.true.)
5060 call readi(controlcard,"NSTEP",n_timestep,1000000)
5061 call readi(controlcard,"NTWE",ntwe,100)
5062 call readi(controlcard,"NTWX",ntwx,1000)
5063 call readi(controlcard,"NFOD",nfodstep,100)
5064 call reada(controlcard,"DT",d_time,1.0d-1)
5065 call reada(controlcard,"DVMAX",dvmax,2.0d1)
5066 call reada(controlcard,"DAMAX",damax,1.0d1)
5067 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
5068 call readi(controlcard,"LANG",lang,0)
5069 RESPA = index(controlcard,"RESPA") .gt. 0
5070 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
5071 ntime_split0=ntime_split
5072 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5073 ntime_split0=ntime_split
5074 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5075 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5076 rest = index(controlcard,"REST").gt.0
5077 tbf = index(controlcard,"TBF").gt.0
5078 usampl = index(controlcard,"USAMPL").gt.0
5079 ! write(iout,*) "KURWA",usampl
5080 mdpdb = index(controlcard,"MDPDB").gt.0
5081 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5082 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5083 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5084 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5085 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5086 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5087 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5088 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5089 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5090 large = index(controlcard,"LARGE").gt.0
5091 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5092 rattle = index(controlcard,"RATTLE").gt.0
5093 preminim=(index(controlcard,'PREMINIM').gt.0)
5094 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5095 write (iout,*) "PREMINIM ",preminim
5096 dccart=(index(controlcard,'CART').gt.0)
5097 if (preminim) call read_minim
5098 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5104 if(me.eq.king.or..not.out1file) then
5106 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5108 write (iout,'(a)') "The units are:"
5109 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5110 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5111 " acceleration: angstrom/(48.9 fs)**2"
5112 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5114 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5115 write (iout,'(a60,f10.5,a)') &
5116 "Initial time step of numerical integration:",d_time,&
5118 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5120 write (iout,'(2a,i4,a)') &
5121 "A-MTS algorithm used; initial time step for fast-varying",&
5122 " short-range forces split into",ntime_split," steps."
5123 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5124 r_cut," lambda",rlamb
5126 write (iout,'(2a,f10.5)') &
5127 "Maximum acceleration threshold to reduce the time step",&
5128 "/increase split number:",damax
5129 write (iout,'(2a,f10.5)') &
5130 "Maximum predicted energy drift to reduce the timestep",&
5131 "/increase split number:",edriftmax
5132 write (iout,'(a60,f10.5)') &
5133 "Maximum velocity threshold to reduce velocities:",dvmax
5134 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5135 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5136 if (rattle) write (iout,'(a60)') &
5137 "Rattle algorithm used to constrain the virtual bonds"
5141 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5142 call reada(controlcard,"RWAT",rwat,1.4d0)
5143 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5144 surfarea=index(controlcard,"SURFAREA").gt.0
5145 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5146 if(me.eq.king.or..not.out1file)then
5147 write (iout,'(/a,$)') "Langevin dynamics calculation"
5149 write (iout,'(a/)') &
5150 " with direct integration of Langevin equations"
5151 else if (lang.eq.2) then
5152 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5153 else if (lang.eq.3) then
5154 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5155 else if (lang.eq.4) then
5156 write (iout,'(a/)') " in overdamped mode"
5158 write (iout,'(//a,i5)') &
5159 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5162 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5163 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5164 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5165 write (iout,'(a60,f10.5)') &
5166 "Scaling factor of the friction forces:",scal_fric
5167 if (surfarea) write (iout,'(2a,i10,a)') &
5168 "Friction coefficients will be scaled by solvent-accessible",&
5169 " surface area every",reset_fricmat," steps."
5171 ! Calculate friction coefficients and bounds of stochastic forces
5172 eta=6*pi*cPoise*etawat
5173 if(me.eq.king.or..not.out1file) &
5174 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5177 do j=1,5 !types of molecules
5178 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5179 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5181 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5182 do j=1,5 !types of molecules
5184 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5185 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5189 if(me.eq.king.or..not.out1file)then
5191 write (iout,'(/2a/)') &
5192 "Radii of site types and friction coefficients and std's of",&
5193 " stochastic forces of fully exposed sites"
5194 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5197 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5198 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5203 if(me.eq.king.or..not.out1file)then
5204 write (iout,'(a)') "Berendsen bath calculation"
5205 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5206 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5208 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5209 count_reset_moment," steps"
5211 write (iout,'(a,i10,a)') &
5212 "Velocities will be reset at random every",count_reset_vel,&
5216 if(me.eq.king.or..not.out1file) &
5217 write (iout,'(a31)') "Microcanonical mode calculation"
5219 if(me.eq.king.or..not.out1file)then
5220 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5222 write(iout,*) "MD running with constraints."
5223 write(iout,*) "Equilibration time ", eq_time, " mtus."
5224 write(iout,*) "Constraining ", nfrag," fragments."
5225 write(iout,*) "Length of each fragment, weight and q0:"
5227 write (iout,*) "Set of restraints #",iset
5229 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5230 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5232 write(iout,*) "constraints between ", npair, "fragments."
5233 write(iout,*) "constraint pairs, weights and q0:"
5235 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5236 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5238 write(iout,*) "angle constraints within ", nfrag_back,&
5239 "backbone fragments."
5240 write(iout,*) "fragment, weights:"
5242 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5243 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5244 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5247 iset=mod(kolor,nset)+1
5250 if(me.eq.king.or..not.out1file) &
5251 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5253 end subroutine read_MDpar
5254 !-----------------------------------------------------------------------------
5258 ! implicit real*8 (a-h,o-z)
5259 ! include 'DIMENSIONS'
5260 ! include 'COMMON.MAP'
5261 ! include 'COMMON.IOUNITS'
5262 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5263 character(len=80) :: mapcard !,ucase
5266 ! real(kind=8) :: var,ene
5269 read (inp,'(a)') mapcard
5270 mapcard=ucase(mapcard)
5271 if (index(mapcard,'PHI').gt.0) then
5273 else if (index(mapcard,'THE').gt.0) then
5275 else if (index(mapcard,'ALP').gt.0) then
5277 else if (index(mapcard,'OME').gt.0) then
5280 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5281 stop 'Error - illegal variable spec in MAP card.'
5283 call readi (mapcard,'RES1',res1(imap),0)
5284 call readi (mapcard,'RES2',res2(imap),0)
5285 if (res1(imap).eq.0) then
5286 res1(imap)=res2(imap)
5287 else if (res2(imap).eq.0) then
5288 res2(imap)=res1(imap)
5290 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5291 write (iout,'(a)') &
5292 'Error - illegal definition of variable group in MAP.'
5293 stop 'Error - illegal definition of variable group in MAP.'
5295 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5296 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5297 call readi(mapcard,'NSTEP',nstep(imap),0)
5298 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5299 write (iout,'(a)') &
5300 'Illegal boundary and/or step size specification in MAP.'
5301 stop 'Illegal boundary and/or step size specification in MAP.'
5305 end subroutine map_read
5306 !-----------------------------------------------------------------------------
5309 use control_data, only: vdisulf
5311 ! implicit real*8 (a-h,o-z)
5312 ! include 'DIMENSIONS'
5313 ! include 'COMMON.IOUNITS'
5314 ! include 'COMMON.GEO'
5315 ! include 'COMMON.CSA'
5316 ! include 'COMMON.BANK'
5317 ! include 'COMMON.CONTROL'
5318 ! character(len=80) :: ucase
5319 character(len=620) :: mcmcard
5321 ! integer :: ntf,ik,iw_pdb
5322 ! real(kind=8) :: var,ene
5324 call card_concat(mcmcard,.true.)
5326 call readi(mcmcard,'NCONF',nconf,50)
5327 call readi(mcmcard,'NADD',nadd,0)
5328 call readi(mcmcard,'JSTART',jstart,1)
5329 call readi(mcmcard,'JEND',jend,1)
5330 call readi(mcmcard,'NSTMAX',nstmax,500000)
5331 call readi(mcmcard,'N0',n0,1)
5332 call readi(mcmcard,'N1',n1,6)
5333 call readi(mcmcard,'N2',n2,4)
5334 call readi(mcmcard,'N3',n3,0)
5335 call readi(mcmcard,'N4',n4,0)
5336 call readi(mcmcard,'N5',n5,0)
5337 call readi(mcmcard,'N6',n6,10)
5338 call readi(mcmcard,'N7',n7,0)
5339 call readi(mcmcard,'N8',n8,0)
5340 call readi(mcmcard,'N9',n9,0)
5341 call readi(mcmcard,'N14',n14,0)
5342 call readi(mcmcard,'N15',n15,0)
5343 call readi(mcmcard,'N16',n16,0)
5344 call readi(mcmcard,'N17',n17,0)
5345 call readi(mcmcard,'N18',n18,0)
5347 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5349 call readi(mcmcard,'NDIFF',ndiff,2)
5350 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5351 call readi(mcmcard,'IS1',is1,1)
5352 call readi(mcmcard,'IS2',is2,8)
5353 call readi(mcmcard,'NRAN0',nran0,4)
5354 call readi(mcmcard,'NRAN1',nran1,2)
5355 call readi(mcmcard,'IRR',irr,1)
5356 call readi(mcmcard,'NSEED',nseed,20)
5357 call readi(mcmcard,'NTOTAL',ntotal,10000)
5358 call reada(mcmcard,'CUT1',cut1,2.0d0)
5359 call reada(mcmcard,'CUT2',cut2,5.0d0)
5360 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5361 call readi(mcmcard,'ICMAX',icmax,3)
5362 call readi(mcmcard,'IRESTART',irestart,0)
5363 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5366 call reada(mcmcard,'DELE',dele,20.0d0)
5367 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5368 call readi(mcmcard,'IREF',iref,0)
5369 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5370 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5371 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5372 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5373 write (iout,*) "NCONF_IN",nconf_in
5375 end subroutine csaread
5376 !-----------------------------------------------------------------------------
5380 use control_data, only: MaxMoveType
5383 ! implicit real*8 (a-h,o-z)
5384 ! include 'DIMENSIONS'
5385 ! include 'COMMON.MCM'
5386 ! include 'COMMON.MCE'
5387 ! include 'COMMON.IOUNITS'
5388 ! character(len=80) :: ucase
5389 character(len=320) :: mcmcard
5392 ! real(kind=8) :: var,ene
5394 call card_concat(mcmcard,.true.)
5395 call readi(mcmcard,'MAXACC',maxacc,100)
5396 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5397 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5398 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5399 call readi(mcmcard,'MAXREPM',maxrepm,200)
5400 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5401 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5402 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5403 call reada(mcmcard,'E_UP',e_up,5.0D0)
5404 call reada(mcmcard,'DELTE',delte,0.1D0)
5405 call readi(mcmcard,'NSWEEP',nsweep,5)
5406 call readi(mcmcard,'NSTEPH',nsteph,0)
5407 call readi(mcmcard,'NSTEPC',nstepc,0)
5408 call reada(mcmcard,'TMIN',tmin,298.0D0)
5409 call reada(mcmcard,'TMAX',tmax,298.0D0)
5410 call readi(mcmcard,'NWINDOW',nwindow,0)
5411 call readi(mcmcard,'PRINT_MC',print_mc,0)
5412 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5413 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5414 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5415 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5416 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5417 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5418 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5419 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5420 if (nwindow.gt.0) then
5421 allocate(winstart(nwindow)) !!el (maxres)
5422 allocate(winend(nwindow)) !!el
5423 allocate(winlen(nwindow)) !!el
5424 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5426 winlen(i)=winend(i)-winstart(i)+1
5429 if (tmax.lt.tmin) tmax=tmin
5430 if (tmax.eq.tmin) then
5434 if (nstepc.gt.0 .and. nsteph.gt.0) then
5435 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5436 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5438 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5439 ! Probabilities of different move types
5440 sumpro_type(0)=0.0D0
5441 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5442 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5443 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5444 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5445 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5446 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5447 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5449 print *,'i',i,' sumprotype',sumpro_type(i)
5450 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5451 print *,'i',i,' sumprotype',sumpro_type(i)
5454 end subroutine mcmread
5455 !-----------------------------------------------------------------------------
5456 subroutine read_minim
5459 ! implicit real*8 (a-h,o-z)
5460 ! include 'DIMENSIONS'
5461 ! include 'COMMON.MINIM'
5462 ! include 'COMMON.IOUNITS'
5463 ! character(len=80) :: ucase
5464 character(len=320) :: minimcard
5466 ! integer :: ntf,ik,iw_pdb
5467 ! real(kind=8) :: var,ene
5469 call card_concat(minimcard,.true.)
5470 call readi(minimcard,'MAXMIN',maxmin,2000)
5471 call readi(minimcard,'MAXFUN',maxfun,5000)
5472 call readi(minimcard,'MINMIN',minmin,maxmin)
5473 call readi(minimcard,'MINFUN',minfun,maxmin)
5474 call reada(minimcard,'TOLF',tolf,1.0D-2)
5475 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5476 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5477 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5478 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5479 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5480 'Options in energy minimization:'
5481 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5482 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5483 'MinMin:',MinMin,' MinFun:',MinFun,&
5484 ' TolF:',TolF,' RTolF:',RTolF
5486 end subroutine read_minim
5487 !-----------------------------------------------------------------------------
5488 subroutine openunits
5490 use MD_data, only: usampl
5493 use control_data, only:out1file
5494 use control, only: getenv_loc
5495 ! implicit real*8 (a-h,o-z)
5496 ! include 'DIMENSIONS'
5499 character(len=16) :: form,nodename
5500 integer :: nodelen,ierror,npos
5502 ! include 'COMMON.SETUP'
5503 ! include 'COMMON.IOUNITS'
5504 ! include 'COMMON.MD'
5505 ! include 'COMMON.CONTROL'
5506 integer :: lenpre,lenpot,lentmp !,ilen
5508 character(len=3) :: out1file_text !,ucase
5509 character(len=3) :: ll
5512 ! integer :: ntf,ik,iw_pdb
5513 ! real(kind=8) :: var,ene
5515 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5516 call getenv_loc("PREFIX",prefix)
5518 call getenv_loc("POT",pot)
5519 call getenv_loc("DIRTMP",tmpdir)
5520 call getenv_loc("CURDIR",curdir)
5521 call getenv_loc("OUT1FILE",out1file_text)
5522 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5523 out1file_text=ucase(out1file_text)
5524 if (out1file_text(1:1).eq."Y") then
5527 out1file=fg_rank.gt.0
5532 if (lentmp.gt.0) then
5533 write (*,'(80(1h!))')
5534 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5535 write (*,'(80(1h!))')
5536 write (*,*)"All output files will be on node /tmp directory."
5538 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5539 if (me.eq.king) then
5540 write (*,*) "The master node is ",nodename
5541 else if (fg_rank.eq.0) then
5542 write (*,*) "I am the CG slave node ",nodename
5544 write (*,*) "I am the FG slave node ",nodename
5547 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5548 lenpre = lentmp+lenpre+1
5550 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5551 ! Get the names and open the input files
5552 #if defined(WINIFL) || defined(WINPGI)
5553 open(1,file=pref_orig(:ilen(pref_orig))// &
5554 '.inp',status='old',readonly,shared)
5555 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5556 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5557 ! Get parameter filenames and open the parameter files.
5558 call getenv_loc('BONDPAR',bondname)
5559 open (ibond,file=bondname,status='old',readonly,shared)
5560 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5561 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5562 call getenv_loc('THETPAR',thetname)
5563 open (ithep,file=thetname,status='old',readonly,shared)
5564 call getenv_loc('ROTPAR',rotname)
5565 open (irotam,file=rotname,status='old',readonly,shared)
5566 call getenv_loc('TORPAR',torname)
5567 open (itorp,file=torname,status='old',readonly,shared)
5568 call getenv_loc('TORDPAR',tordname)
5569 open (itordp,file=tordname,status='old',readonly,shared)
5570 call getenv_loc('FOURIER',fouriername)
5571 open (ifourier,file=fouriername,status='old',readonly,shared)
5572 call getenv_loc('ELEPAR',elename)
5573 open (ielep,file=elename,status='old',readonly,shared)
5574 call getenv_loc('SIDEPAR',sidename)
5575 open (isidep,file=sidename,status='old',readonly,shared)
5577 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5578 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5579 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5580 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5581 call getenv_loc('TORPAR_NUCL',torname_nucl)
5582 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5583 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5584 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5585 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5586 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5589 #elif (defined CRAY) || (defined AIX)
5590 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5592 ! print *,"Processor",myrank," opened file 1"
5593 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5594 ! print *,"Processor",myrank," opened file 9"
5595 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5596 ! Get parameter filenames and open the parameter files.
5597 call getenv_loc('BONDPAR',bondname)
5598 open (ibond,file=bondname,status='old',action='read')
5599 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5600 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5602 ! print *,"Processor",myrank," opened file IBOND"
5603 call getenv_loc('THETPAR',thetname)
5604 open (ithep,file=thetname,status='old',action='read')
5605 ! print *,"Processor",myrank," opened file ITHEP"
5606 call getenv_loc('ROTPAR',rotname)
5607 open (irotam,file=rotname,status='old',action='read')
5608 ! print *,"Processor",myrank," opened file IROTAM"
5609 call getenv_loc('TORPAR',torname)
5610 open (itorp,file=torname,status='old',action='read')
5611 ! print *,"Processor",myrank," opened file ITORP"
5612 call getenv_loc('TORDPAR',tordname)
5613 open (itordp,file=tordname,status='old',action='read')
5614 ! print *,"Processor",myrank," opened file ITORDP"
5615 call getenv_loc('SCCORPAR',sccorname)
5616 open (isccor,file=sccorname,status='old',action='read')
5617 ! print *,"Processor",myrank," opened file ISCCOR"
5618 call getenv_loc('FOURIER',fouriername)
5619 open (ifourier,file=fouriername,status='old',action='read')
5620 ! print *,"Processor",myrank," opened file IFOURIER"
5621 call getenv_loc('ELEPAR',elename)
5622 open (ielep,file=elename,status='old',action='read')
5623 ! print *,"Processor",myrank," opened file IELEP"
5624 call getenv_loc('SIDEPAR',sidename)
5625 open (isidep,file=sidename,status='old',action='read')
5627 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5628 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5629 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5630 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5631 call getenv_loc('TORPAR_NUCL',torname_nucl)
5632 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5633 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5634 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5635 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5636 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5638 call getenv_loc('LIPTRANPAR',liptranname)
5639 open (iliptranpar,file=liptranname,status='old',action='read')
5640 call getenv_loc('TUBEPAR',tubename)
5641 open (itube,file=tubename,status='old',action='read')
5642 call getenv_loc('IONPAR',ionname)
5643 open (iion,file=ionname,status='old',action='read')
5644 call getenv_loc('IONPAR_TRAN',iontranname)
5645 open (iiontran,file=iontranname,status='old',action='read')
5646 ! print *,"Processor",myrank," opened file ISIDEP"
5647 ! print *,"Processor",myrank," opened parameter files"
5649 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5650 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5651 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5652 ! Get parameter filenames and open the parameter files.
5653 call getenv_loc('BONDPAR',bondname)
5654 open (ibond,file=bondname,status='old')
5655 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5656 open (ibond_nucl,file=bondname_nucl,status='old')
5658 call getenv_loc('THETPAR',thetname)
5659 open (ithep,file=thetname,status='old')
5660 call getenv_loc('ROTPAR',rotname)
5661 open (irotam,file=rotname,status='old')
5662 call getenv_loc('TORPAR',torname)
5663 open (itorp,file=torname,status='old')
5664 call getenv_loc('TORDPAR',tordname)
5665 open (itordp,file=tordname,status='old')
5666 call getenv_loc('SCCORPAR',sccorname)
5667 open (isccor,file=sccorname,status='old')
5668 call getenv_loc('FOURIER',fouriername)
5669 open (ifourier,file=fouriername,status='old')
5670 call getenv_loc('ELEPAR',elename)
5671 open (ielep,file=elename,status='old')
5672 call getenv_loc('SIDEPAR',sidename)
5673 open (isidep,file=sidename,status='old')
5675 open (ithep_nucl,file=thetname_nucl,status='old')
5676 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5677 open (irotam_nucl,file=rotname_nucl,status='old')
5678 call getenv_loc('TORPAR_NUCL',torname_nucl)
5679 open (itorp_nucl,file=torname_nucl,status='old')
5680 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5681 open (itordp_nucl,file=tordname_nucl,status='old')
5682 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5683 open (isidep_nucl,file=sidename_nucl,status='old')
5685 call getenv_loc('LIPTRANPAR',liptranname)
5686 open (iliptranpar,file=liptranname,status='old')
5687 call getenv_loc('TUBEPAR',tubename)
5688 open (itube,file=tubename,status='old')
5689 call getenv_loc('IONPAR',ionname)
5690 open (iion,file=ionname,status='old')
5691 call getenv_loc('IONPAR_NUCL',ionnuclname)
5692 open (iionnucl,file=ionnuclname,status='old')
5693 call getenv_loc('IONPAR_TRAN',iontranname)
5694 open (iiontran,file=iontranname,status='old')
5695 call getenv_loc('WATWAT',iwaterwatername)
5696 open (iwaterwater,file=iwaterwatername,status='old')
5697 call getenv_loc('WATPROT',iwaterscname)
5698 open (iwatersc,file=iwaterscname,status='old')
5701 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5703 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5704 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5705 ! Get parameter filenames and open the parameter files.
5706 call getenv_loc('BONDPAR',bondname)
5707 open (ibond,file=bondname,status='old',action='read')
5708 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5709 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5710 call getenv_loc('THETPAR',thetname)
5711 open (ithep,file=thetname,status='old',action='read')
5712 call getenv_loc('ROTPAR',rotname)
5713 open (irotam,file=rotname,status='old',action='read')
5714 call getenv_loc('TORPAR',torname)
5715 open (itorp,file=torname,status='old',action='read')
5716 call getenv_loc('TORDPAR',tordname)
5717 open (itordp,file=tordname,status='old',action='read')
5718 call getenv_loc('SCCORPAR',sccorname)
5719 open (isccor,file=sccorname,status='old',action='read')
5721 call getenv_loc('THETPARPDB',thetname_pdb)
5722 print *,"thetname_pdb ",thetname_pdb
5723 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5724 print *,ithep_pdb," opened"
5726 call getenv_loc('FOURIER',fouriername)
5727 open (ifourier,file=fouriername,status='old',readonly)
5728 call getenv_loc('ELEPAR',elename)
5729 open (ielep,file=elename,status='old',readonly)
5730 call getenv_loc('SIDEPAR',sidename)
5731 open (isidep,file=sidename,status='old',readonly)
5733 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5734 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5735 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5736 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5737 call getenv_loc('TORPAR_NUCL',torname_nucl)
5738 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5739 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5740 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5741 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5742 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5743 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5744 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5745 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5746 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5747 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5748 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5749 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5750 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5753 call getenv_loc('LIPTRANPAR',liptranname)
5754 open (iliptranpar,file=liptranname,status='old',action='read')
5755 call getenv_loc('LIPBOND',lipbondname)
5756 open (ilipbond,file=lipbondname,status='old',action='read')
5757 call getenv_loc('LIPNONBOND',lipnonbondname)
5758 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5759 call getenv_loc('TUBEPAR',tubename)
5760 open (itube,file=tubename,status='old',action='read')
5761 call getenv_loc('IONPAR',ionname)
5762 open (iion,file=ionname,status='old',action='read')
5763 call getenv_loc('IONPAR_NUCL',ionnuclname)
5764 open (iionnucl,file=ionnuclname,status='old',action='read')
5765 call getenv_loc('IONPAR_TRAN',iontranname)
5766 open (iiontran,file=iontranname,status='old',action='read')
5767 call getenv_loc('WATWAT',iwaterwatername)
5768 open (iwaterwater,file=iwaterwatername,status='old',action='read')
5769 call getenv_loc('WATPROT',iwaterscname)
5770 open (iwatersc,file=iwaterscname,status='old',action='read')
5773 call getenv_loc('ROTPARPDB',rotname_pdb)
5774 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5777 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5778 #if defined(WINIFL) || defined(WINPGI)
5779 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5780 #elif (defined CRAY) || (defined AIX)
5781 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5783 open (iscpp_nucl,file=scpname_nucl,status='old')
5785 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5790 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5791 ! Use -DOLDSCP to use hard-coded constants instead.
5793 call getenv_loc('SCPPAR',scpname)
5794 #if defined(WINIFL) || defined(WINPGI)
5795 open (iscpp,file=scpname,status='old',readonly,shared)
5796 #elif (defined CRAY) || (defined AIX)
5797 open (iscpp,file=scpname,status='old',action='read')
5799 open (iscpp,file=scpname,status='old')
5801 open (iscpp,file=scpname,status='old',action='read')
5804 call getenv_loc('PATTERN',patname)
5805 #if defined(WINIFL) || defined(WINPGI)
5806 open (icbase,file=patname,status='old',readonly,shared)
5807 #elif (defined CRAY) || (defined AIX)
5808 open (icbase,file=patname,status='old',action='read')
5810 open (icbase,file=patname,status='old')
5812 open (icbase,file=patname,status='old',action='read')
5815 ! Open output file only for CG processes
5816 ! print *,"Processor",myrank," fg_rank",fg_rank
5817 if (fg_rank.eq.0) then
5819 if (nodes.eq.1) then
5822 npos = dlog10(dfloat(nodes-1))+1
5824 if (npos.lt.3) npos=3
5825 write (liczba,'(i1)') npos
5826 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5828 write (liczba,form) me
5829 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5830 liczba(:ilen(liczba))
5831 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5833 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5835 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5836 liczba(:ilen(liczba))//'.mol2'
5837 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5838 liczba(:ilen(liczba))//'.stat'
5840 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5841 //liczba(:ilen(liczba))//'.stat')
5842 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5845 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5846 liczba(:ilen(liczba))//'.const'
5851 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5852 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5853 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5854 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5855 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5857 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5859 rest2name=prefix(:ilen(prefix))//'.rst'
5861 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5864 #if defined(AIX) || defined(PGI)
5865 if (me.eq.king .or. .not. out1file) &
5866 open(iout,file=outname,status='unknown')
5868 if (fg_rank.gt.0) then
5869 write (liczba,'(i3.3)') myrank/nfgtasks
5870 write (ll,'(bz,i3.3)') fg_rank
5871 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5876 open(igeom,file=intname,status='unknown',position='append')
5877 open(ipdb,file=pdbname,status='unknown')
5878 open(imol2,file=mol2name,status='unknown')
5879 open(istat,file=statname,status='unknown',position='append')
5881 !1out open(iout,file=outname,status='unknown')
5884 if (me.eq.king .or. .not.out1file) &
5885 open(iout,file=outname,status='unknown')
5887 if (fg_rank.gt.0) then
5888 write (liczba,'(i3.3)') myrank/nfgtasks
5889 write (ll,'(bz,i3.3)') fg_rank
5890 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5895 open(igeom,file=intname,status='unknown',access='append')
5896 open(ipdb,file=pdbname,status='unknown')
5897 open(imol2,file=mol2name,status='unknown')
5898 open(istat,file=statname,status='unknown',access='append')
5900 !1out open(iout,file=outname,status='unknown')
5903 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5904 csa_seed=prefix(:lenpre)//'.CSA.seed'
5905 csa_history=prefix(:lenpre)//'.CSA.history'
5906 csa_bank=prefix(:lenpre)//'.CSA.bank'
5907 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5908 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5909 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5910 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5911 csa_int=prefix(:lenpre)//'.int'
5912 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5913 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5914 csa_in=prefix(:lenpre)//'.CSA.in'
5915 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5918 write (iout,'(80(1h-))')
5919 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5920 write (iout,'(80(1h-))')
5921 write (iout,*) "Input file : ",&
5922 pref_orig(:ilen(pref_orig))//'.inp'
5923 write (iout,*) "Output file : ",&
5924 outname(:ilen(outname))
5926 write (iout,*) "Sidechain potential file : ",&
5927 sidename(:ilen(sidename))
5929 write (iout,*) "SCp potential file : ",&
5930 scpname(:ilen(scpname))
5932 write (iout,*) "Electrostatic potential file : ",&
5933 elename(:ilen(elename))
5934 write (iout,*) "Cumulant coefficient file : ",&
5935 fouriername(:ilen(fouriername))
5936 write (iout,*) "Torsional parameter file : ",&
5937 torname(:ilen(torname))
5938 write (iout,*) "Double torsional parameter file : ",&
5939 tordname(:ilen(tordname))
5940 write (iout,*) "SCCOR parameter file : ",&
5941 sccorname(:ilen(sccorname))
5942 write (iout,*) "Bond & inertia constant file : ",&
5943 bondname(:ilen(bondname))
5944 write (iout,*) "Bending parameter file : ",&
5945 thetname(:ilen(thetname))
5946 write (iout,*) "Rotamer parameter file : ",&
5947 rotname(:ilen(rotname))
5950 write (iout,*) "Thetpdb parameter file : ",&
5951 thetname_pdb(:ilen(thetname_pdb))
5954 write (iout,*) "Threading database : ",&
5955 patname(:ilen(patname))
5957 write (iout,*)" DIRTMP : ",&
5959 write (iout,'(80(1h-))')
5962 end subroutine openunits
5963 !-----------------------------------------------------------------------------
5966 use geometry_data, only: nres,dc
5968 ! implicit real*8 (a-h,o-z)
5969 ! include 'DIMENSIONS'
5970 ! include 'COMMON.CHAIN'
5971 ! include 'COMMON.IOUNITS'
5972 ! include 'COMMON.MD'
5975 ! real(kind=8) :: var,ene
5977 open(irest2,file=rest2name,status='unknown')
5978 read(irest2,*) totT,EK,potE,totE,t_bath
5981 ! AL 4/17/17: Now reading d_t(0,:) too
5983 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5986 ! AL 4/17/17: Now reading d_c(0,:) too
5988 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5991 read (irest2,*) iset
5995 end subroutine readrst
5996 !-----------------------------------------------------------------------------
5997 subroutine read_fragments
6001 use control_data, only:out1file
6004 ! implicit real*8 (a-h,o-z)
6005 ! include 'DIMENSIONS'
6009 ! include 'COMMON.SETUP'
6010 ! include 'COMMON.CHAIN'
6011 ! include 'COMMON.IOUNITS'
6012 ! include 'COMMON.MD'
6013 ! include 'COMMON.CONTROL'
6016 ! real(kind=8) :: var,ene
6018 read(inp,*) nset,nfrag,npair,nfrag_back
6020 !el from module energy
6021 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
6022 if(.not.allocated(wfrag_back)) then
6023 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6024 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6026 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
6027 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
6029 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
6030 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
6033 if(me.eq.king.or..not.out1file) &
6034 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
6035 " nfrag_back",nfrag_back
6037 read(inp,*) mset(iset)
6039 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
6041 if(me.eq.king.or..not.out1file) &
6042 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
6043 ifrag(2,i,iset), qinfrag(i,iset)
6046 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
6048 if(me.eq.king.or..not.out1file) &
6049 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
6050 ipair(2,i,iset), qinpair(i,iset)
6053 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6054 wfrag_back(3,i,iset),&
6055 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6056 if(me.eq.king.or..not.out1file) &
6057 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6058 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6062 end subroutine read_fragments
6063 !-----------------------------------------------------------------------------
6065 !-----------------------------------------------------------------------------
6069 ! implicit real*8 (a-h,o-z)
6070 ! include 'DIMENSIONS'
6071 ! include 'COMMON.CSA'
6072 ! include 'COMMON.BANK'
6073 ! include 'COMMON.IOUNITS'
6075 ! integer :: ntf,ik,iw_pdb
6076 ! real(kind=8) :: var,ene
6078 open(icsa_in,file=csa_in,status="old",err=100)
6079 read(icsa_in,*) nconf
6080 read(icsa_in,*) jstart,jend
6081 read(icsa_in,*) nstmax
6082 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6083 read(icsa_in,*) nran0,nran1,irr
6084 read(icsa_in,*) nseed
6085 read(icsa_in,*) ntotal,cut1,cut2
6086 read(icsa_in,*) estop
6087 read(icsa_in,*) icmax,irestart
6088 read(icsa_in,*) ntbankm,dele,difcut
6089 read(icsa_in,*) iref,rmscut,pnccut
6090 read(icsa_in,*) ndiff
6097 end subroutine csa_read
6098 !-----------------------------------------------------------------------------
6099 subroutine initial_write
6102 ! implicit real*8 (a-h,o-z)
6103 ! include 'DIMENSIONS'
6104 ! include 'COMMON.CSA'
6105 ! include 'COMMON.BANK'
6106 ! include 'COMMON.IOUNITS'
6108 ! integer :: ntf,ik,iw_pdb
6109 ! real(kind=8) :: var,ene
6111 open(icsa_seed,file=csa_seed,status="unknown")
6112 write(icsa_seed,*) "seed"
6114 #if defined(AIX) || defined(PGI)
6115 open(icsa_history,file=csa_history,status="unknown",&
6118 open(icsa_history,file=csa_history,status="unknown",&
6121 write(icsa_history,*) nconf
6122 write(icsa_history,*) jstart,jend
6123 write(icsa_history,*) nstmax
6124 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6125 write(icsa_history,*) nran0,nran1,irr
6126 write(icsa_history,*) nseed
6127 write(icsa_history,*) ntotal,cut1,cut2
6128 write(icsa_history,*) estop
6129 write(icsa_history,*) icmax,irestart
6130 write(icsa_history,*) ntbankm,dele,difcut
6131 write(icsa_history,*) iref,rmscut,pnccut
6132 write(icsa_history,*) ndiff
6134 write(icsa_history,*)
6137 open(icsa_bank1,file=csa_bank1,status="unknown")
6138 write(icsa_bank1,*) 0
6142 end subroutine initial_write
6143 !-----------------------------------------------------------------------------
6144 subroutine restart_write
6147 ! implicit real*8 (a-h,o-z)
6148 ! include 'DIMENSIONS'
6149 ! include 'COMMON.IOUNITS'
6150 ! include 'COMMON.CSA'
6151 ! include 'COMMON.BANK'
6153 ! integer :: ntf,ik,iw_pdb
6154 ! real(kind=8) :: var,ene
6156 #if defined(AIX) || defined(PGI)
6157 open(icsa_history,file=csa_history,position="append")
6159 open(icsa_history,file=csa_history,access="append")
6161 write(icsa_history,*)
6162 write(icsa_history,*) "This is restart"
6163 write(icsa_history,*)
6164 write(icsa_history,*) nconf
6165 write(icsa_history,*) jstart,jend
6166 write(icsa_history,*) nstmax
6167 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6168 write(icsa_history,*) nran0,nran1,irr
6169 write(icsa_history,*) nseed
6170 write(icsa_history,*) ntotal,cut1,cut2
6171 write(icsa_history,*) estop
6172 write(icsa_history,*) icmax,irestart
6173 write(icsa_history,*) ntbankm,dele,difcut
6174 write(icsa_history,*) iref,rmscut,pnccut
6175 write(icsa_history,*) ndiff
6176 write(icsa_history,*)
6177 write(icsa_history,*) "irestart is: ", irestart
6179 write(icsa_history,*)
6183 end subroutine restart_write
6184 !-----------------------------------------------------------------------------
6186 !-----------------------------------------------------------------------------
6187 subroutine write_pdb(npdb,titelloc,ee)
6189 ! implicit real*8 (a-h,o-z)
6190 ! include 'DIMENSIONS'
6191 ! include 'COMMON.IOUNITS'
6192 character(len=50) :: titelloc1
6193 character*(*) :: titelloc
6194 character(len=3) :: zahl
6195 character(len=5) :: liczba5
6197 integer :: npdb !,ilen
6201 ! real(kind=8) :: var,ene
6205 if (npdb.lt.1000) then
6206 call numstr(npdb,zahl)
6207 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6209 if (npdb.lt.10000) then
6210 write(liczba5,'(i1,i4)') 0,npdb
6212 write(liczba5,'(i5)') npdb
6214 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6216 call pdbout(ee,titelloc1,ipdb)
6219 end subroutine write_pdb
6220 !-----------------------------------------------------------------------------
6222 !-----------------------------------------------------------------------------
6223 subroutine write_thread_summary
6224 ! Thread the sequence through a database of known structures
6225 use control_data, only: refstr
6227 use energy_data, only: n_ene_comp
6229 ! implicit real*8 (a-h,o-z)
6230 ! include 'DIMENSIONS'
6232 use MPI_data !include 'COMMON.INFO'
6235 ! include 'COMMON.CONTROL'
6236 ! include 'COMMON.CHAIN'
6237 ! include 'COMMON.DBASE'
6238 ! include 'COMMON.INTERACT'
6239 ! include 'COMMON.VAR'
6240 ! include 'COMMON.THREAD'
6241 ! include 'COMMON.FFIELD'
6242 ! include 'COMMON.SBRIDGE'
6243 ! include 'COMMON.HEADER'
6244 ! include 'COMMON.NAMES'
6245 ! include 'COMMON.IOUNITS'
6246 ! include 'COMMON.TIME1'
6248 integer,dimension(maxthread) :: ip
6249 real(kind=8),dimension(0:n_ene) :: energia
6251 integer :: i,j,ii,jj,ipj,ik,kk,ist
6252 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6254 write (iout,'(30x,a/)') &
6255 ' *********** Summary threading statistics ************'
6256 write (iout,'(a)') 'Initial energies:'
6257 write (iout,'(a4,2x,a12,14a14,3a8)') &
6258 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6259 'RMSnat','NatCONT','NNCONT','RMS'
6260 ! Energy sort patterns
6265 enet=ener(n_ene-1,ip(i))
6268 if (ener(n_ene-1,ip(j)).lt.enet) then
6270 enet=ener(n_ene-1,ip(j))
6282 ist=nres_base(2,ii)+ipatt(2,i)
6284 energia(i)=ener0(kk,i)
6286 etot=ener0(n_ene_comp+1,i)
6287 rmsnat=ener0(n_ene_comp+2,i)
6288 rms=ener0(n_ene_comp+3,i)
6289 frac=ener0(n_ene_comp+4,i)
6290 frac_nn=ener0(n_ene_comp+5,i)
6293 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6294 i,str_nam(ii),ist+1,&
6295 (energia(print_order(kk)),kk=1,nprint_ene),&
6296 etot,rmsnat,frac,frac_nn,rms
6298 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6299 i,str_nam(ii),ist+1,&
6300 (energia(print_order(kk)),kk=1,nprint_ene),etot
6303 write (iout,'(//a)') 'Final energies:'
6304 write (iout,'(a4,2x,a12,17a14,3a8)') &
6305 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6306 'RMSnat','NatCONT','NNCONT','RMS'
6310 ist=nres_base(2,ii)+ipatt(2,i)
6312 energia(kk)=ener(kk,ik)
6314 etot=ener(n_ene_comp+1,i)
6315 rmsnat=ener(n_ene_comp+2,i)
6316 rms=ener(n_ene_comp+3,i)
6317 frac=ener(n_ene_comp+4,i)
6318 frac_nn=ener(n_ene_comp+5,i)
6319 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6320 i,str_nam(ii),ist+1,&
6321 (energia(print_order(kk)),kk=1,nprint_ene),&
6322 etot,rmsnat,frac,frac_nn,rms
6324 write (iout,'(/a/)') 'IEXAM array:'
6325 write (iout,'(i5)') nexcl
6327 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6329 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6330 'Max. time for threading step ',max_time_for_thread,&
6331 'Average time for threading step: ',ave_time_for_thread
6333 end subroutine write_thread_summary
6334 !-----------------------------------------------------------------------------
6335 subroutine write_stat_thread(ithread,ipattern,ist)
6337 use energy_data, only: n_ene_comp
6339 ! implicit real*8 (a-h,o-z)
6340 ! include "DIMENSIONS"
6341 ! include "COMMON.CONTROL"
6342 ! include "COMMON.IOUNITS"
6343 ! include "COMMON.THREAD"
6344 ! include "COMMON.FFIELD"
6345 ! include "COMMON.DBASE"
6346 ! include "COMMON.NAMES"
6347 real(kind=8),dimension(0:n_ene) :: energia
6349 integer :: ithread,ipattern,ist,i
6350 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6352 #if defined(AIX) || defined(PGI)
6353 open(istat,file=statname,position='append')
6355 open(istat,file=statname,access='append')
6358 energia(i)=ener(i,ithread)
6360 etot=ener(n_ene_comp+1,ithread)
6361 rmsnat=ener(n_ene_comp+2,ithread)
6362 rms=ener(n_ene_comp+3,ithread)
6363 frac=ener(n_ene_comp+4,ithread)
6364 frac_nn=ener(n_ene_comp+5,ithread)
6365 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6366 ithread,str_nam(ipattern),ist+1,&
6367 (energia(print_order(i)),i=1,nprint_ene),&
6368 etot,rmsnat,frac,frac_nn,rms
6371 end subroutine write_stat_thread
6372 !-----------------------------------------------------------------------------
6374 subroutine readpdb_template(k)
6375 ! Read the PDB file for read_constr_homology with read2sigma
6376 ! and convert the peptide geometry into virtual-chain geometry.
6378 ! include 'DIMENSIONS'
6379 ! include 'COMMON.LOCAL'
6380 ! include 'COMMON.VAR'
6381 ! include 'COMMON.CHAIN'
6382 ! include 'COMMON.INTERACT'
6383 ! include 'COMMON.IOUNITS'
6384 ! include 'COMMON.GEO'
6385 ! include 'COMMON.NAMES'
6386 ! include 'COMMON.CONTROL'
6387 ! include 'COMMON.FRAG'
6388 ! include 'COMMON.SETUP'
6389 use compare_data, only:nhfrag,nbfrag
6390 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6392 logical lprn /.false./,fail
6393 real(kind=8), dimension (3):: e1,e2,e3
6394 real(kind=8) :: dcj,efree_temp
6398 real(kind=8), dimension (3,20) :: sccor
6400 integer, dimension (:), allocatable :: iterter
6401 if(.not.allocated(iterter))allocate(iterter(nres))
6408 write (2,*) "UNRES_PDB",unres_pdb
6416 read (ipdbin,'(a80)',end=10) card
6417 if (card(:3).eq.'END') then
6419 else if (card(:3).eq.'TER') then
6422 itype(ires_old-1,1)=ntyp1
6423 iterter(ires_old-1)=1
6424 #if defined(WHAM_RUN) || defined(CLUSTER)
6425 if (ires_old.lt.nres) then
6427 itype(ires_old,1)=ntyp1
6429 #if defined(WHAM_RUN) || defined(CLUSTER)
6433 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6436 dc(j,ires)=sccor(j,iii)
6439 call sccenter(ires,iii,sccor)
6442 ! Fish out the ATOM cards.
6443 if (index(card(1:4),'ATOM').gt.0) then
6444 read (card(12:16),*) atom
6445 ! write (iout,*) "! ",atom," !",ires
6446 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6447 read (card(23:26),*) ires
6448 read (card(18:20),'(a3)') res
6449 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6450 ! & " ires_old",ires_old
6451 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6452 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6453 if (ires-ishift+ishift1.ne.ires_old) then
6454 ! Calculate the CM of the preceding residue.
6458 dc(j,ires_old)=sccor(j,iii)
6461 call sccenter(ires_old,iii,sccor)
6465 ! Start new residue.
6466 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6469 else if (ibeg.eq.1) then
6470 ! write (iout,*) "BEG ires",ires
6472 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6476 ires=ires-ishift+ishift1
6478 ! write (iout,*) "ishift",ishift," ires",ires,
6479 ! & " ires_old",ires_old
6480 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6482 else if (ibeg.eq.2) then
6484 ishift=-ires_old+ires-1
6486 ! write (iout,*) "New chain started",ires,ishift
6489 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6490 ires=ires-ishift+ishift1
6493 if (res.eq.'ACE' .or. res.eq.'NHE') then
6496 itype(ires,1)=rescode(ires,res,0,1)
6499 ires=ires-ishift+ishift1
6501 ! write (iout,*) "ires_old",ires_old," ires",ires
6502 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6505 ! write (2,*) "ires",ires," res ",res," ity",ity
6506 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6507 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6508 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6509 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6511 write (iout,'(2i3,2x,a,3f8.3)') &
6512 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6516 sccor(j,iii)=c(j,ires)
6518 if (ishift.ne.0) then
6519 ires_ca=ires+ishift-ishift1
6523 ! write (*,*) card(23:27),ires,itype(ires)
6524 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6525 atom.ne.'N' .and. atom.ne.'C' .and.&
6526 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6527 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6528 ! write (iout,*) "sidechain ",atom
6530 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6534 10 if(me.eq.king.or..not.out1file) &
6535 write (iout,'(a,i5)') ' Nres: ',ires
6536 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6539 write(2,*) "tutaj",ires,nres
6541 ! write (iout,*) i,itype(i),itype(i+1)
6542 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6543 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6544 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6545 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6546 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6548 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6549 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6556 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6560 dcj=(c(j,i-2)-c(j,i-3))/2.0
6561 if (dcj.eq.0) dcj=1.23591524223
6566 else !itype(i+1).eq.ntyp1
6568 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6569 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6576 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6580 dcj=(c(j,i+3)-c(j,i+2))/2.0
6581 if (dcj.eq.0) dcj=1.23591524223
6586 endif !itype(i+1).eq.ntyp1
6587 endif !itype.eq.ntyp1
6589 ! Calculate the CM of the last side chain.
6592 dc(j,ires)=sccor(j,iii)
6595 call sccenter(ires,iii,sccor)
6599 if (itype(nres,1).ne.10) then
6603 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6604 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6611 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6615 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6616 if (dcj.eq.0) dcj=1.23591524223
6617 c(j,nres)=c(j,nres-1)+dcj
6618 c(j,2*nres)=c(j,nres)
6629 c(j,2*nres)=c(j,nres)
6631 if (itype(1,1).eq.ntyp1) then
6635 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6636 call refsys(2,3,4,e1,e2,e3,fail)
6643 c(j,1)=c(j,2)-1.9d0*e2(j)
6647 dcj=(c(j,4)-c(j,3))/2.0
6653 ! Copy the coordinates to reference coordinates
6659 ! Calculate internal coordinates.
6660 if (out_template_coord) then
6661 write (iout,'(/a)') &
6662 "Cartesian coordinates of the reference structure"
6663 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6664 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6666 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6667 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6668 (c(j,ires+nres),j=1,3)
6671 ! Calculate internal coordinates.
6672 call int_from_cart(.true.,out_template_coord)
6673 call sc_loc_geom(.false.)
6675 thetaref(i)=theta(i)
6680 dc(j,i)=c(j,i+1)-c(j,i)
6681 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6686 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6687 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6689 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6690 ! & vbld_inv(i+nres)
6695 cref(j,i+nres,1)=c(j,i+nres)
6705 end subroutine readpdb_template
6706 !-----------------------------------------------------------------------------
6708 !-----------------------------------------------------------------------------
6709 end module io_config