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
534 write (iout,*) 'FTORS factor =',ftors
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
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii,ncatprotparm
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
780 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
781 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(nbondterm(ntyp)) !(ntyp)
783 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
784 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(msc(ntyp+1,5)) !(ntyp+1)
786 allocate(isc(ntyp+1,5)) !(ntyp+1)
787 allocate(restok(ntyp+1,5)) !(ntyp+1)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 read (ibond,*) vbldp0,akp,mp,ip,pstok
798 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
799 dsc(i) = vbldsc0(1,i)
803 dsc_inv(i)=1.0D0/dsc(i)
807 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
809 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
810 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
811 dsc(i) = vbldsc0(1,i)
815 dsc_inv(i)=1.0D0/dsc(i)
819 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
822 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
823 ! dsc(i) = vbldsc0_nucl(1,i)
827 ! dsc_inv(i)=1.0D0/dsc(i)
830 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
831 ! do i=1,ntyp_molec(2)
832 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
833 ! aksc_nucl(j,i),abond0_nucl(j,i),&
834 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
835 ! dsc(i) = vbldsc0(1,i)
839 ! dsc_inv(i)=1.0D0/dsc(i)
844 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
845 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
847 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
849 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
850 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
852 write (iout,'(13x,3f10.5)') &
853 vbldsc0(j,i),aksc(j,i),abond0(j,i)
858 read(iion,*) msc(i,5),restok(i,5)
859 print *,msc(i,5),restok(i,5)
863 read (iion,*) ncatprotparm
864 allocate(catprm(ncatprotparm,4))
866 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
869 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
870 !----------------------------------------------------
871 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
872 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
873 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
874 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
875 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
876 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
886 allocate(liptranene(ntyp))
887 !C reading lipid parameters
888 write (iout,*) "iliptranpar",iliptranpar
890 read(iliptranpar,*) pepliptran
893 read(iliptranpar,*) liptranene(i)
894 print *,liptranene(i)
900 ! Read the parameters of the probability distribution/energy expression
901 ! of the virtual-bond valence angles theta
904 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
905 (bthet(j,i,1,1),j=1,2)
906 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
907 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
908 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
912 athet(1,i,1,-1)=athet(1,i,1,1)
913 athet(2,i,1,-1)=athet(2,i,1,1)
914 bthet(1,i,1,-1)=-bthet(1,i,1,1)
915 bthet(2,i,1,-1)=-bthet(2,i,1,1)
916 athet(1,i,-1,1)=-athet(1,i,1,1)
917 athet(2,i,-1,1)=-athet(2,i,1,1)
918 bthet(1,i,-1,1)=bthet(1,i,1,1)
919 bthet(2,i,-1,1)=bthet(2,i,1,1)
923 athet(1,i,-1,-1)=athet(1,-i,1,1)
924 athet(2,i,-1,-1)=-athet(2,-i,1,1)
925 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
926 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
927 athet(1,i,-1,1)=athet(1,-i,1,1)
928 athet(2,i,-1,1)=-athet(2,-i,1,1)
929 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
930 bthet(2,i,-1,1)=bthet(2,-i,1,1)
931 athet(1,i,1,-1)=-athet(1,-i,1,1)
932 athet(2,i,1,-1)=athet(2,-i,1,1)
933 bthet(1,i,1,-1)=bthet(1,-i,1,1)
934 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
939 polthet(j,i)=polthet(j,-i)
942 gthet(j,i)=gthet(j,-i)
950 'Parameters of the virtual-bond valence angles:'
951 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
952 ' ATHETA0 ',' A1 ',' A2 ',&
955 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
956 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
958 write (iout,'(/a/9x,5a/79(1h-))') &
959 'Parameters of the expression for sigma(theta_c):',&
960 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
961 ' ALPH3 ',' SIGMA0C '
963 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
964 (polthet(j,i),j=0,3),sigc0(i)
966 write (iout,'(/a/9x,5a/79(1h-))') &
967 'Parameters of the second gaussian:',&
968 ' THETA0 ',' SIGMA0 ',' G1 ',&
971 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
972 sig0(i),(gthet(j,i),j=1,3)
976 'Parameters of the virtual-bond valence angles:'
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Coefficients of expansion',&
979 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
980 ' b1*10^1 ',' b2*10^1 '
982 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
983 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
984 (10*bthet(j,i,1,1),j=1,2)
986 write (iout,'(/a/9x,5a/79(1h-))') &
987 'Parameters of the expression for sigma(theta_c):',&
988 ' alpha0 ',' alph1 ',' alph2 ',&
989 ' alhp3 ',' sigma0c '
991 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
992 (polthet(j,i),j=0,3),sigc0(i)
994 write (iout,'(/a/9x,5a/79(1h-))') &
995 'Parameters of the second gaussian:',&
996 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
999 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1000 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1006 ! Read the parameters of Utheta determined from ab initio surfaces
1007 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1009 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1010 ntheterm3,nsingle,ndouble
1011 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1013 !----------------------------------------------------
1014 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1015 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1016 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1017 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1018 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1019 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1020 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1021 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1022 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1023 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1024 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1025 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1026 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1029 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1030 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1031 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1035 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1037 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1039 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1041 ithetyp(i)=-ithetyp(-i)
1044 aa0thet(:,:,:,:)=0.0d0
1045 aathet(:,:,:,:,:)=0.0d0
1046 bbthet(:,:,:,:,:,:)=0.0d0
1047 ccthet(:,:,:,:,:,:)=0.0d0
1048 ddthet(:,:,:,:,:,:)=0.0d0
1049 eethet(:,:,:,:,:,:)=0.0d0
1050 ffthet(:,:,:,:,:,:,:)=0.0d0
1051 ggthet(:,:,:,:,:,:,:)=0.0d0
1053 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1055 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1056 ! VAR:1=non-glicyne non-proline 2=proline
1057 ! VAR:negative values for D-aminoacid
1059 do j=-nthetyp,nthetyp
1060 do k=-nthetyp,nthetyp
1061 read (ithep,'(6a)',end=111,err=111) res1
1062 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1063 ! VAR: aa0thet is variable describing the average value of Foureir
1064 ! VAR: expansion series
1065 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1066 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1067 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1068 read (ithep,*,end=111,err=111) &
1069 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1070 read (ithep,*,end=111,err=111) &
1071 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1072 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1073 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1074 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1076 read (ithep,*,end=111,err=111) &
1077 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1078 ffthet(lll,llll,ll,i,j,k,iblock),&
1079 ggthet(llll,lll,ll,i,j,k,iblock),&
1080 ggthet(lll,llll,ll,i,j,k,iblock),&
1081 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1086 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1087 ! coefficients of theta-and-gamma-dependent terms are zero.
1088 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1089 ! RECOMENTDED AFTER VERSION 3.3)
1093 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1094 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1096 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1097 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1100 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1102 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1105 ! AND COMMENT THE LOOPS BELOW
1109 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1110 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1112 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1113 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1116 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1118 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1123 ! Substitution for D aminoacids from symmetry.
1126 do j=-nthetyp,nthetyp
1127 do k=-nthetyp,nthetyp
1128 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1130 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1134 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1135 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1136 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1137 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1143 ffthet(llll,lll,ll,i,j,k,iblock)= &
1144 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1145 ffthet(lll,llll,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(llll,lll,ll,-i,-j,-k,iblock)
1149 ggthet(lll,llll,ll,i,j,k,iblock)= &
1150 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1159 ! Control printout of the coefficients of virtual-bond-angle potentials
1162 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1167 write (iout,'(//4a)') &
1168 'Type ',onelett(i),onelett(j),onelett(k)
1169 write (iout,'(//a,10x,a)') " l","a[l]"
1170 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1171 write (iout,'(i2,1pe15.5)') &
1172 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1174 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1175 "b",l,"c",l,"d",l,"e",l
1177 write (iout,'(i2,4(1pe15.5))') m,&
1178 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1179 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1183 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1184 "f+",l,"f-",l,"g+",l,"g-",l
1187 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1188 ffthet(n,m,l,i,j,k,iblock),&
1189 ffthet(m,n,l,i,j,k,iblock),&
1190 ggthet(n,m,l,i,j,k,iblock),&
1191 ggthet(m,n,l,i,j,k,iblock)
1201 write (2,*) "Start reading THETA_PDB",ithep_pdb
1203 ! write (2,*) 'i=',i
1204 read (ithep_pdb,*,err=111,end=111) &
1205 a0thet(i),(athet(j,i,1,1),j=1,2),&
1206 (bthet(j,i,1,1),j=1,2)
1207 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1208 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1209 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1210 sigc0(i)=sigc0(i)**2
1213 athet(1,i,1,-1)=athet(1,i,1,1)
1214 athet(2,i,1,-1)=athet(2,i,1,1)
1215 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1216 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1217 athet(1,i,-1,1)=-athet(1,i,1,1)
1218 athet(2,i,-1,1)=-athet(2,i,1,1)
1219 bthet(1,i,-1,1)=bthet(1,i,1,1)
1220 bthet(2,i,-1,1)=bthet(2,i,1,1)
1223 a0thet(i)=a0thet(-i)
1224 athet(1,i,-1,-1)=athet(1,-i,1,1)
1225 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1226 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1227 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1228 athet(1,i,-1,1)=athet(1,-i,1,1)
1229 athet(2,i,-1,1)=-athet(2,-i,1,1)
1230 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1232 athet(1,i,1,-1)=-athet(1,-i,1,1)
1233 athet(2,i,1,-1)=athet(2,-i,1,1)
1234 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1235 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1236 theta0(i)=theta0(-i)
1240 polthet(j,i)=polthet(j,-i)
1243 gthet(j,i)=gthet(j,-i)
1246 write (2,*) "End reading THETA_PDB"
1250 !--------------- Reading theta parameters for nucleic acid-------
1251 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1252 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1253 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1254 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1255 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1256 nthetyp_nucl+1,nthetyp_nucl+1))
1257 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1258 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1259 nthetyp_nucl+1,nthetyp_nucl+1))
1260 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1261 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1262 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1263 nthetyp_nucl+1,nthetyp_nucl+1))
1264 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1265 nthetyp_nucl+1,nthetyp_nucl+1))
1266 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1269 nthetyp_nucl+1,nthetyp_nucl+1))
1270 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1271 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1272 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1273 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1274 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1275 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1277 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1278 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1280 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1282 aa0thet_nucl(:,:,:)=0.0d0
1283 aathet_nucl(:,:,:,:)=0.0d0
1284 bbthet_nucl(:,:,:,:,:)=0.0d0
1285 ccthet_nucl(:,:,:,:,:)=0.0d0
1286 ddthet_nucl(:,:,:,:,:)=0.0d0
1287 eethet_nucl(:,:,:,:,:)=0.0d0
1288 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1289 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1294 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1295 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1296 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1297 read (ithep_nucl,*,end=111,err=111) &
1298 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1299 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1300 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1301 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1302 read (ithep_nucl,*,end=111,err=111) &
1303 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1304 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1305 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1310 !-------------------------------------------
1311 allocate(nlob(ntyp1)) !(ntyp1)
1312 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1313 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1314 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1321 gaussc(:,:,:,:)=0.0D0
1325 ! Read the parameters of the probability distribution/energy expression
1326 ! of the side chains.
1329 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1333 dsc_inv(i)=1.0D0/dsc(i)
1344 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1345 ((blower(k,l,1),l=1,k),k=1,3)
1346 censc(1,1,-i)=censc(1,1,i)
1347 censc(2,1,-i)=censc(2,1,i)
1348 censc(3,1,-i)=-censc(3,1,i)
1350 read (irotam,*,end=112,err=112) bsc(j,i)
1351 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1352 ((blower(k,l,j),l=1,k),k=1,3)
1353 censc(1,j,-i)=censc(1,j,i)
1354 censc(2,j,-i)=censc(2,j,i)
1355 censc(3,j,-i)=-censc(3,j,i)
1356 ! BSC is amplitude of Gaussian
1363 akl=akl+blower(k,m,j)*blower(l,m,j)
1367 if (((k.eq.3).and.(l.ne.3)) &
1368 .or.((l.eq.3).and.(k.ne.3))) then
1369 gaussc(k,l,j,-i)=-akl
1370 gaussc(l,k,j,-i)=-akl
1372 gaussc(k,l,j,-i)=akl
1373 gaussc(l,k,j,-i)=akl
1382 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1385 if (nlobi.gt.0) then
1387 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1388 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1389 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1390 'log h',(bsc(j,i),j=1,nlobi)
1391 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1392 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1394 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1395 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1398 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1399 write (iout,'(a,f10.4,4(16x,f10.4))') &
1400 'Center ',(bsc(j,i),j=1,nlobi)
1401 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1410 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1411 ! added by Urszula Kozlowska 07/11/2007
1413 !el Maximum number of SC local term fitting function coefficiants
1414 !el integer,parameter :: maxsccoef=65
1416 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1419 read (irotam,*,end=112,err=112)
1421 read (irotam,*,end=112,err=112)
1424 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1428 !---------reading nucleic acid parameters for rotamers-------------------
1429 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1430 do i=1,ntyp_molec(2)
1431 read (irotam_nucl,*,end=112,err=112)
1433 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1439 write (iout,*) "Base rotamer parameters"
1440 do i=1,ntyp_molec(2)
1441 write (iout,'(a)') restyp(i,2)
1442 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1447 ! Read the parameters of the probability distribution/energy expression
1448 ! of the side chains.
1450 write (2,*) "Start reading ROTAM_PDB"
1452 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1456 dsc_inv(i)=1.0D0/dsc(i)
1467 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1468 ((blower(k,l,1),l=1,k),k=1,3)
1470 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1471 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1472 ((blower(k,l,j),l=1,k),k=1,3)
1479 akl=akl+blower(k,m,j)*blower(l,m,j)
1489 write (2,*) "End reading ROTAM_PDB"
1495 ! Read torsional parameters in old format
1497 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1499 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1500 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1501 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1503 !el from energy module--------
1504 allocate(v1(nterm_old,ntortyp,ntortyp))
1505 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1506 !el---------------------------
1511 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1517 write (iout,'(/a/)') 'Torsional constants:'
1520 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1521 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1527 ! Read torsional parameters
1529 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1530 read (itorp,*,end=113,err=113) ntortyp
1531 !el from energy module---------
1532 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1533 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1535 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1536 allocate(vlor2(maxlor,ntortyp,ntortyp))
1537 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1538 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1540 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1541 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1542 !el---------------------------
1545 !el---------------------------
1547 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1549 itortyp(i)=-itortyp(-i)
1551 itortyp(ntyp1)=ntortyp
1552 itortyp(-ntyp1)=-ntortyp
1554 write (iout,*) 'ntortyp',ntortyp
1556 do j=-ntortyp+1,ntortyp-1
1557 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1559 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1560 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1563 do k=1,nterm(i,j,iblock)
1564 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1566 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1567 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1568 v0ij=v0ij+si*v1(k,i,j,iblock)
1570 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1571 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1572 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1574 do k=1,nlor(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1576 vlor2(k,i,j),vlor3(k,i,j)
1577 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1580 v0(-i,-j,iblock)=v0ij
1586 write (iout,'(/a/)') 'Torsional constants:'
1588 do i=-ntortyp,ntortyp
1589 do j=-ntortyp,ntortyp
1590 write (iout,*) 'ityp',i,' jtyp',j
1591 write (iout,*) 'Fourier constants'
1592 do k=1,nterm(i,j,iblock)
1593 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1596 write (iout,*) 'Lorenz constants'
1597 do k=1,nlor(i,j,iblock)
1598 write (iout,'(3(1pe15.5))') &
1599 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1605 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1607 ! 6/23/01 Read parameters for double torsionals
1609 !el from energy module------------
1610 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1611 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1612 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1613 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1614 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1615 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1616 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1617 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1618 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1619 !---------------------------------
1623 do j=-ntortyp+1,ntortyp-1
1624 do k=-ntortyp+1,ntortyp-1
1625 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1626 ! write (iout,*) "OK onelett",
1629 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1630 .or. t3.ne.toronelet(k)) then
1631 write (iout,*) "Error in double torsional parameter file",&
1634 call MPI_Finalize(Ierror)
1636 stop "Error in double torsional parameter file"
1638 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1639 ntermd_2(i,j,k,iblock)
1640 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1641 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1642 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1643 ntermd_1(i,j,k,iblock))
1644 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1645 ntermd_1(i,j,k,iblock))
1646 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1647 ntermd_1(i,j,k,iblock))
1648 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1649 ntermd_1(i,j,k,iblock))
1650 ! Martix of D parameters for one dimesional foureir series
1651 do l=1,ntermd_1(i,j,k,iblock)
1652 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1653 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1654 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1655 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1656 ! write(iout,*) "whcodze" ,
1657 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1659 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1660 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1661 v2s(m,l,i,j,k,iblock),&
1662 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1663 ! Martix of D parameters for two dimesional fourier series
1664 do l=1,ntermd_2(i,j,k,iblock)
1666 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1667 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1668 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1669 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1678 write (iout,*) 'Constants for double torsionals'
1681 do j=-ntortyp+1,ntortyp-1
1682 do k=-ntortyp+1,ntortyp-1
1683 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1684 ' nsingle',ntermd_1(i,j,k,iblock),&
1685 ' ndouble',ntermd_2(i,j,k,iblock)
1687 write (iout,*) 'Single angles:'
1688 do l=1,ntermd_1(i,j,k,iblock)
1689 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1690 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1691 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1692 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1695 write (iout,*) 'Pairs of angles:'
1696 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1697 do l=1,ntermd_2(i,j,k,iblock)
1698 write (iout,'(i5,20f10.5)') &
1699 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1702 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1703 do l=1,ntermd_2(i,j,k,iblock)
1704 write (iout,'(i5,20f10.5)') &
1705 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1706 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1715 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1716 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1717 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1718 !el from energy module---------
1719 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1720 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1722 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1723 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1724 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1725 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1727 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1728 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1729 !el---------------------------
1732 !el--------------------
1733 read (itorp_nucl,*,end=113,err=113) &
1734 (itortyp_nucl(i),i=1,ntyp_molec(2))
1735 ! print *,itortyp_nucl(:)
1736 !c write (iout,*) 'ntortyp',ntortyp
1739 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1740 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1743 do k=1,nterm_nucl(i,j)
1744 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1745 v0ij=v0ij+si*v1_nucl(k,i,j)
1748 do k=1,nlor_nucl(i,j)
1749 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1750 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1751 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1757 ! Read of Side-chain backbone correlation parameters
1758 ! Modified 11 May 2012 by Adasko
1761 read (isccor,*,end=119,err=119) nsccortyp
1763 !el from module energy-------------
1764 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1765 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1766 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1767 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1768 !-----------------------------------
1770 !el from module energy-------------
1771 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1773 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1775 isccortyp(i)=-isccortyp(-i)
1777 iscprol=isccortyp(20)
1778 ! write (iout,*) 'ntortyp',ntortyp
1780 !c maxinter is maximum interaction sites
1781 !el from module energy---------
1782 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1783 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1784 -nsccortyp:nsccortyp))
1785 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1786 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1787 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1788 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1789 !-----------------------------------
1793 read (isccor,*,end=119,err=119) &
1794 nterm_sccor(i,j),nlor_sccor(i,j)
1800 nterm_sccor(-i,j)=nterm_sccor(i,j)
1801 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1802 nterm_sccor(i,-j)=nterm_sccor(i,j)
1803 do k=1,nterm_sccor(i,j)
1804 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1806 if (j.eq.iscprol) then
1807 if (i.eq.isccortyp(10)) then
1808 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1809 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1812 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1813 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1814 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1815 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1816 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1817 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1818 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1821 if (i.eq.isccortyp(10)) then
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1823 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1825 if (j.eq.isccortyp(10)) then
1826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1829 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1830 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1831 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1832 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1833 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1834 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1838 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1839 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1840 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1841 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1844 do k=1,nlor_sccor(i,j)
1845 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1846 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1847 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1848 (1+vlor3sccor(k,i,j)**2)
1850 v0sccor(l,i,j)=v0ijsccor
1851 v0sccor(l,-i,j)=v0ijsccor1
1852 v0sccor(l,i,-j)=v0ijsccor2
1853 v0sccor(l,-i,-j)=v0ijsccor3
1859 !el from module energy-------------
1860 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1862 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1863 ! write (iout,*) 'ntortyp',ntortyp
1865 !c maxinter is maximum interaction sites
1866 !el from module energy---------
1867 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1868 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1869 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1870 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1871 !-----------------------------------
1875 read (isccor,*,end=119,err=119) &
1876 nterm_sccor(i,j),nlor_sccor(i,j)
1880 do k=1,nterm_sccor(i,j)
1881 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1883 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1886 do k=1,nlor_sccor(i,j)
1887 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1888 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1889 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1890 (1+vlor3sccor(k,i,j)**2)
1892 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1900 write (iout,'(/a/)') 'Torsional constants:'
1903 write (iout,*) 'ityp',i,' jtyp',j
1904 write (iout,*) 'Fourier constants'
1905 do k=1,nterm_sccor(i,j)
1906 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1908 write (iout,*) 'Lorenz constants'
1909 do k=1,nlor_sccor(i,j)
1910 write (iout,'(3(1pe15.5))') &
1911 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1918 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1919 ! interaction energy of the Gly, Ala, and Pro prototypes.
1924 write (iout,*) "Coefficients of the cumulants"
1926 read (ifourier,*) nloctyp
1927 !write(iout,*) "nloctyp",nloctyp
1928 !el from module energy-------
1929 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1930 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1931 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1932 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1933 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1934 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1935 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1936 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1940 !--------------------------------
1943 read (ifourier,*,end=115,err=115)
1944 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1946 write (iout,*) 'Type',i
1947 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1957 B1tilde(1,-i) =-b(3)
1959 ! b1tilde(1,i)=0.0d0
1960 ! b1tilde(2,i)=0.0d0
1985 Ctilde(1,2,-i)=-b(9)
1989 ! Ctilde(1,1,i)=0.0d0
1990 ! Ctilde(1,2,i)=0.0d0
1991 ! Ctilde(2,1,i)=0.0d0
1992 ! Ctilde(2,2,i)=0.0d0
2010 Dtilde(1,2,-i)=-b(8)
2014 ! Dtilde(1,1,i)=0.0d0
2015 ! Dtilde(1,2,i)=0.0d0
2016 ! Dtilde(2,1,i)=0.0d0
2017 ! Dtilde(2,2,i)=0.0d0
2018 EE(1,1,i)= b(10)+b(11)
2019 EE(2,2,i)=-b(10)+b(11)
2020 EE(2,1,i)= b(12)-b(13)
2021 EE(1,2,i)= b(12)+b(13)
2022 EE(1,1,-i)= b(10)+b(11)
2023 EE(2,2,-i)=-b(10)+b(11)
2024 EE(2,1,-i)=-b(12)+b(13)
2025 EE(1,2,-i)=-b(12)-b(13)
2031 ! ee(2,1,i)=ee(1,2,i)
2035 write (iout,*) 'Type',i
2037 write(iout,*) B1(1,i),B1(2,i)
2039 write(iout,*) B2(1,i),B2(2,i)
2042 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2046 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2050 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2055 ! Read electrostatic-interaction parameters
2060 write (iout,'(/a)') 'Electrostatic interaction constants:'
2061 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2062 'IT','JT','APP','BPP','AEL6','AEL3'
2064 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2065 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2066 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2067 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2072 app (i,j)=epp(i,j)*rri*rri
2073 bpp (i,j)=-2.0D0*epp(i,j)*rri
2074 ael6(i,j)=elpp6(i,j)*4.2D0**6
2075 ael3(i,j)=elpp3(i,j)*4.2D0**3
2077 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2083 ! Read side-chain interaction parameters.
2085 !el from module energy - COMMON.INTERACT-------
2086 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2087 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2088 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2090 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2091 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2092 allocate(epslip(ntyp,ntyp))
2100 !--------------------------------
2102 read (isidep,*,end=117,err=117) ipot,expon
2103 if (ipot.lt.1 .or. ipot.gt.5) then
2104 write (iout,'(2a)') 'Error while reading SC interaction',&
2105 'potential file - unknown potential type.'
2107 call MPI_Finalize(Ierror)
2113 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2114 ', exponents are ',expon,2*expon
2115 ! goto (10,20,30,30,40) ipot
2117 !----------------------- LJ potential ---------------------------------
2119 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2120 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121 (sigma0(i),i=1,ntyp)
2123 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2124 write (iout,'(a/)') 'The epsilon array:'
2125 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2126 write (iout,'(/a)') 'One-body parameters:'
2127 write (iout,'(a,4x,a)') 'residue','sigma'
2128 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2131 !----------------------- LJK potential --------------------------------
2133 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2134 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2135 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2137 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2138 write (iout,'(a/)') 'The epsilon array:'
2139 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2140 write (iout,'(/a)') 'One-body parameters:'
2141 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2142 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2146 !---------------------- GB or BP potential -----------------------------
2150 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2152 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2153 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2154 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2155 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2157 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2160 ! For the GB potential convert sigma'**2 into chi'
2163 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2167 write (iout,'(/a/)') 'Parameters of the BP potential:'
2168 write (iout,'(a/)') 'The epsilon array:'
2169 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2170 write (iout,'(/a)') 'One-body parameters:'
2171 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2173 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2174 chip(i),alp(i),i=1,ntyp)
2177 !--------------------- GBV potential -----------------------------------
2179 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2180 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2181 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2182 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2184 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2185 write (iout,'(a/)') 'The epsilon array:'
2186 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2187 write (iout,'(/a)') 'One-body parameters:'
2188 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2189 's||/s_|_^2',' chip ',' alph '
2190 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2191 sigii(i),chip(i),alp(i),i=1,ntyp)
2194 write(iout,*)"Wrong ipot"
2200 !-----------------------------------------------------------------------
2201 ! Calculate the "working" parameters of SC interactions.
2203 !el from module energy - COMMON.INTERACT-------
2204 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2205 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2206 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2207 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2209 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2222 sc_aa_tube_par(:)=0.0d0
2223 sc_bb_tube_par(:)=0.0d0
2225 !--------------------------------
2230 epslip(i,j)=epslip(j,i)
2235 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2236 sigma(j,i)=sigma(i,j)
2237 rs0(i,j)=dwa16*sigma(i,j)
2241 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2242 'Working parameters of the SC interactions:',&
2243 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2248 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2257 sigeps=dsign(1.0D0,epsij)
2259 aa_aq(i,j)=epsij*rrij*rrij
2260 bb_aq(i,j)=-sigeps*epsij*rrij
2261 aa_aq(j,i)=aa_aq(i,j)
2262 bb_aq(j,i)=bb_aq(i,j)
2263 epsijlip=epslip(i,j)
2264 sigeps=dsign(1.0D0,epsijlip)
2265 epsijlip=dabs(epsijlip)
2266 aa_lip(i,j)=epsijlip*rrij*rrij
2267 bb_lip(i,j)=-sigeps*epsijlip*rrij
2268 aa_lip(j,i)=aa_lip(i,j)
2269 bb_lip(j,i)=bb_lip(i,j)
2270 !C write(iout,*) aa_lip
2272 sigt1sq=sigma0(i)**2
2273 sigt2sq=sigma0(j)**2
2276 ratsig1=sigt2sq/sigt1sq
2277 ratsig2=1.0D0/ratsig1
2278 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2279 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2280 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2284 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2285 sigmaii(i,j)=rsum_max
2286 sigmaii(j,i)=rsum_max
2288 ! sigmaii(i,j)=r0(i,j)
2289 ! sigmaii(j,i)=r0(i,j)
2291 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2292 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2293 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2294 augm(i,j)=epsij*r_augm**(2*expon)
2295 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2302 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2303 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2304 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2309 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2310 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2311 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2312 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2313 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2314 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2315 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2316 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2317 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2318 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2320 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2322 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2324 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2333 read (isidep_nucl,*) ipot_nucl
2334 ! print *,"TU?!",ipot_nucl
2335 if (ipot_nucl.eq.1) then
2336 do i=1,ntyp_molec(2)
2337 do j=i,ntyp_molec(2)
2338 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2339 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2343 do i=1,ntyp_molec(2)
2344 do j=i,ntyp_molec(2)
2345 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2346 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2347 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2351 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2352 do i=1,ntyp_molec(2)
2353 do j=i,ntyp_molec(2)
2354 rrij=sigma_nucl(i,j)
2358 epsij=4*eps_nucl(i,j)
2359 sigeps=dsign(1.0D0,epsij)
2361 aa_nucl(i,j)=epsij*rrij*rrij
2362 bb_nucl(i,j)=-sigeps*epsij*rrij
2363 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2364 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2365 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2366 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2367 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2368 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2371 aa_nucl(i,j)=aa_nucl(j,i)
2372 bb_nucl(i,j)=bb_nucl(j,i)
2373 ael3_nucl(i,j)=ael3_nucl(j,i)
2374 ael6_nucl(i,j)=ael6_nucl(j,i)
2375 ael63_nucl(i,j)=ael63_nucl(j,i)
2376 ael32_nucl(i,j)=ael32_nucl(j,i)
2377 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2378 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2379 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2380 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2381 eps_nucl(i,j)=eps_nucl(j,i)
2382 sigma_nucl(i,j)=sigma_nucl(j,i)
2383 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2387 write(iout,*) "tube param"
2388 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2389 ccavtubpep,dcavtubpep,tubetranenepep
2390 sigmapeptube=sigmapeptube**6
2391 sigeps=dsign(1.0D0,epspeptube)
2392 epspeptube=dabs(epspeptube)
2393 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2394 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2395 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2397 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2398 ccavtub(i),dcavtub(i),tubetranene(i)
2399 sigmasctube=sigmasctube**6
2400 sigeps=dsign(1.0D0,epssctube)
2401 epssctube=dabs(epssctube)
2402 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2403 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2404 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2406 !-----------------READING SC BASE POTENTIALS-----------------------------
2407 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2408 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2409 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2410 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2411 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2412 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2413 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2414 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2415 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2416 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2417 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2418 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2419 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2420 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2421 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2422 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2425 do i=1,ntyp_molec(1)
2426 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2427 write (*,*) "Im in ", i, " ", j
2428 read(isidep_scbase,*) &
2429 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2430 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2431 write(*,*) "eps",eps_scbase(i,j)
2432 read(isidep_scbase,*) &
2433 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2434 chis_scbase(i,j,1),chis_scbase(i,j,2)
2435 read(isidep_scbase,*) &
2436 dhead_scbasei(i,j), &
2437 dhead_scbasej(i,j), &
2438 rborn_scbasei(i,j),rborn_scbasej(j,i)
2439 read(isidep_scbase,*) &
2440 (wdipdip_scbase(k,i,j),k=1,3), &
2441 (wqdip_scbase(k,i,j),k=1,2)
2442 read(isidep_scbase,*) &
2443 alphapol_scbase(i,j), &
2444 epsintab_scbase(i,j)
2447 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2448 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2450 do i=1,ntyp_molec(1)
2451 do j=1,ntyp_molec(2)-1
2452 epsij=eps_scbase(i,j)
2453 rrij=sigma_scbase(i,j)
2458 sigeps=dsign(1.0D0,epsij)
2460 aa_scbase(i,j)=epsij*rrij*rrij
2461 bb_scbase(i,j)=-sigeps*epsij*rrij
2464 !-----------------READING PEP BASE POTENTIALS-------------------
2465 allocate(eps_pepbase(ntyp_molec(2)))
2466 allocate(sigma_pepbase(ntyp_molec(2)))
2467 allocate(chi_pepbase(ntyp_molec(2),2))
2468 allocate(chipp_pepbase(ntyp_molec(2),2))
2469 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2470 allocate(sigmap1_pepbase(ntyp_molec(2)))
2471 allocate(sigmap2_pepbase(ntyp_molec(2)))
2472 allocate(chis_pepbase(ntyp_molec(2),2))
2473 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2476 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2477 write (*,*) "Im in ", i, " ", j
2478 read(isidep_pepbase,*) &
2479 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2480 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2481 write(*,*) "eps",eps_pepbase(j)
2482 read(isidep_pepbase,*) &
2483 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2484 chis_pepbase(j,1),chis_pepbase(j,2)
2485 read(isidep_pepbase,*) &
2486 (wdipdip_pepbase(k,j),k=1,3)
2488 allocate(aa_pepbase(ntyp_molec(2)))
2489 allocate(bb_pepbase(ntyp_molec(2)))
2491 do j=1,ntyp_molec(2)-1
2492 epsij=eps_pepbase(j)
2493 rrij=sigma_pepbase(j)
2498 sigeps=dsign(1.0D0,epsij)
2500 aa_pepbase(j)=epsij*rrij*rrij
2501 bb_pepbase(j)=-sigeps*epsij*rrij
2503 !--------------READING SC PHOSPHATE-------------------------------------
2504 allocate(eps_scpho(ntyp_molec(1)))
2505 allocate(sigma_scpho(ntyp_molec(1)))
2506 allocate(chi_scpho(ntyp_molec(1),2))
2507 allocate(chipp_scpho(ntyp_molec(1),2))
2508 allocate(alphasur_scpho(4,ntyp_molec(1)))
2509 allocate(sigmap1_scpho(ntyp_molec(1)))
2510 allocate(sigmap2_scpho(ntyp_molec(1)))
2511 allocate(chis_scpho(ntyp_molec(1),2))
2512 allocate(wqq_scpho(ntyp_molec(1)))
2513 allocate(wqdip_scpho(2,ntyp_molec(1)))
2514 allocate(alphapol_scpho(ntyp_molec(1)))
2515 allocate(epsintab_scpho(ntyp_molec(1)))
2516 allocate(dhead_scphoi(ntyp_molec(1)))
2517 allocate(rborn_scphoi(ntyp_molec(1)))
2518 allocate(rborn_scphoj(ntyp_molec(1)))
2521 do j=1,ntyp_molec(1) ! without U then we will take T for U
2522 write (*,*) "Im in scpho ", i, " ", j
2523 read(isidep_scpho,*) &
2524 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2525 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2526 write(*,*) "eps",eps_scpho(j)
2527 read(isidep_scpho,*) &
2528 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2529 chis_scpho(j,1),chis_scpho(j,2)
2530 read(isidep_scpho,*) &
2531 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2532 read(isidep_scpho,*) &
2533 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j)
2536 allocate(aa_scpho(ntyp_molec(1)))
2537 allocate(bb_scpho(ntyp_molec(1)))
2539 do j=1,ntyp_molec(1)
2546 sigeps=dsign(1.0D0,epsij)
2548 aa_scpho(j)=epsij*rrij*rrij
2549 bb_scpho(j)=-sigeps*epsij*rrij
2553 read(isidep_peppho,*) &
2554 eps_peppho,sigma_peppho
2555 read(isidep_peppho,*) &
2556 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2557 read(isidep_peppho,*) &
2558 (wqdip_peppho(k),k=1,2)
2566 sigeps=dsign(1.0D0,epsij)
2568 aa_peppho=epsij*rrij*rrij
2569 bb_peppho=-sigeps*epsij*rrij
2572 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2577 ! Define the SC-p interaction constants (hard-coded; old style)
2580 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2582 ! aad(i,1)=0.3D0*4.0D0**12
2583 ! Following line for constants currently implemented
2584 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2585 aad(i,1)=1.5D0*4.0D0**12
2586 ! aad(i,1)=0.17D0*5.6D0**12
2588 ! "Soft" SC-p repulsion
2590 ! Following line for constants currently implemented
2591 ! aad(i,1)=0.3D0*4.0D0**6
2592 ! "Hard" SC-p repulsion
2593 bad(i,1)=3.0D0*4.0D0**6
2594 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2603 ! 8/9/01 Read the SC-p interaction constants from file
2606 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2609 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2610 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2611 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2612 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2616 write (iout,*) "Parameters of SC-p interactions:"
2618 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2619 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2624 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2626 do i=1,ntyp_molec(2)
2627 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2629 do i=1,ntyp_molec(2)
2630 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2631 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2633 r0pp=1.12246204830937298142*5.16158
2639 ! Define the constants of the disulfide bridge
2643 ! Old arbitrary potential - commented out.
2648 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2649 ! energy surface of diethyl disulfide.
2650 ! A. Liwo and U. Kozlowska, 11/24/03
2667 write (iout,'(/a)') "Disulfide bridge parameters:"
2668 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2669 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2670 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2671 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2675 111 write (iout,*) "Error reading bending energy parameters."
2677 112 write (iout,*) "Error reading rotamer energy parameters."
2679 113 write (iout,*) "Error reading torsional energy parameters."
2681 114 write (iout,*) "Error reading double torsional energy parameters."
2683 115 write (iout,*) &
2684 "Error reading cumulant (multibody energy) parameters."
2686 116 write (iout,*) "Error reading electrostatic energy parameters."
2688 117 write (iout,*) "Error reading side chain interaction parameters."
2690 118 write (iout,*) "Error reading SCp interaction parameters."
2692 119 write (iout,*) "Error reading SCCOR parameters"
2695 call MPI_Finalize(Ierror)
2699 end subroutine parmread
2701 !-----------------------------------------------------------------------------
2703 !-----------------------------------------------------------------------------
2704 subroutine printmat(ldim,m,n,iout,key,a)
2707 character(len=3),dimension(n) :: key
2708 real(kind=8),dimension(ldim,n) :: a
2710 integer :: i,j,k,m,iout,nlim
2714 write (iout,1000) (key(k),k=i,nlim)
2716 1000 format (/5x,8(6x,a3))
2717 1020 format (/80(1h-)/)
2719 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2722 1010 format (a3,2x,8(f9.4))
2724 end subroutine printmat
2725 !-----------------------------------------------------------------------------
2727 !-----------------------------------------------------------------------------
2729 ! Read the PDB file and convert the peptide geometry into virtual-chain
2732 use energy_data, only: itype,istype
2736 use control, only: rescode,sugarcode
2737 ! implicit real*8 (a-h,o-z)
2738 ! include 'DIMENSIONS'
2739 ! include 'COMMON.LOCAL'
2740 ! include 'COMMON.VAR'
2741 ! include 'COMMON.CHAIN'
2742 ! include 'COMMON.INTERACT'
2743 ! include 'COMMON.IOUNITS'
2744 ! include 'COMMON.GEO'
2745 ! include 'COMMON.NAMES'
2746 ! include 'COMMON.CONTROL'
2747 ! include 'COMMON.DISTFIT'
2748 ! include 'COMMON.SETUP'
2749 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2751 logical :: lprn=.true.,fail
2752 real(kind=8),dimension(3) :: e1,e2,e3
2753 real(kind=8) :: dcj,efree_temp
2754 character(len=3) :: seq,res
2755 character(len=5) :: atom
2756 character(len=80) :: card
2757 real(kind=8),dimension(3,20) :: sccor
2758 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2759 integer :: isugar,molecprev,firstion
2760 character*1 :: sugar
2762 real(kind=8),dimension(3) :: ccc
2764 integer,dimension(2,maxres/3) :: hfrag_alloc
2765 integer,dimension(4,maxres/3) :: bfrag_alloc
2766 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2767 real(kind=8),dimension(:,:), allocatable :: c_temporary
2768 integer,dimension(:,:) , allocatable :: itype_temporary
2769 integer,dimension(:),allocatable :: istype_temp
2776 ! write (2,*) "UNRES_PDB",unres_pdb
2784 !-----------------------------
2785 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2786 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2789 read (ipdbin,'(a80)',end=10) card
2790 ! write (iout,'(a)') card
2791 if (card(:5).eq.'HELIX') then
2794 read(card(22:25),*) hfrag(1,nhfrag)
2795 read(card(34:37),*) hfrag(2,nhfrag)
2797 if (card(:5).eq.'SHEET') then
2800 read(card(24:26),*) bfrag(1,nbfrag)
2801 read(card(35:37),*) bfrag(2,nbfrag)
2802 !rc----------------------------------------
2803 !rc to be corrected !!!
2804 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2805 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2806 !rc----------------------------------------
2808 if (card(:3).eq.'END') then
2810 else if (card(:3).eq.'TER') then
2815 itype(ires_old,molecule)=ntyp1_molec(molecule)
2816 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2817 nres_molec(molecule)=nres_molec(molecule)+2
2819 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2822 dc(j,ires)=sccor(j,iii)
2825 call sccenter(ires,iii,sccor)
2831 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2832 ! Fish out the ATOM cards.
2833 ! write(iout,*) 'card',card(1:20)
2834 if (index(card(1:4),'ATOM').gt.0) then
2835 read (card(12:16),*) atom
2836 ! write (iout,*) "! ",atom," !",ires
2837 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2838 read (card(23:26),*) ires
2839 read (card(18:20),'(a3)') res
2840 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2841 ! & " ires_old",ires_old
2842 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2843 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2844 if (ires-ishift+ishift1.ne.ires_old) then
2845 ! Calculate the CM of the preceding residue.
2846 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2848 ! write (iout,*) "Calculating sidechain center iii",iii
2851 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2854 call sccenter(ires_old,iii,sccor)
2858 ! Start new residue.
2859 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2862 else if (ibeg.eq.1) then
2863 write (iout,*) "BEG ires",ires
2865 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2868 nres_molec(molecule)=nres_molec(molecule)+1
2870 ires=ires-ishift+ishift1
2872 ! write (iout,*) "ishift",ishift," ires",ires,&
2873 ! " ires_old",ires_old
2875 else if (ibeg.eq.2) then
2877 ishift=-ires_old+ires-1 !!!!!
2878 ishift1=ishift1-1 !!!!!
2879 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2880 ires=ires-ishift+ishift1
2881 ! print *,ires,ishift,ishift1
2885 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2886 ires=ires-ishift+ishift1
2889 ! print *,'atom',ires,atom
2890 if (res.eq.'ACE' .or. res.eq.'NHE') then
2893 if (atom.eq.'CA '.or.atom.eq.'N ') then
2895 itype(ires,molecule)=rescode(ires,res,0,molecule)
2897 ! nres_molec(molecule)=nres_molec(molecule)+1
2900 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2901 ! nres_molec(molecule)=nres_molec(molecule)+1
2902 read (card(19:19),'(a1)') sugar
2903 isugar=sugarcode(sugar,ires)
2904 ! if (ibeg.eq.1) then
2908 ! print *,"ires=",ires,istype(ires)
2914 ires=ires-ishift+ishift1
2916 ! write (iout,*) "ires_old",ires_old," ires",ires
2917 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2920 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2921 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2922 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2923 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2924 ! print *,ires,ishift,ishift1
2925 ! write (iout,*) "backbone ",atom
2927 write (iout,'(2i3,2x,a,3f8.3)') &
2928 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2931 nres_molec(molecule)=nres_molec(molecule)+1
2933 sccor(j,iii)=c(j,ires)
2935 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2936 atom.eq."C2'" .or. atom.eq."C3'" &
2937 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2938 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2939 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2940 ! print *,ires,ishift,ishift1
2944 ! sccor(j,iii)=c(j,ires)
2947 c(j,ires)=c(j,ires)+ccc(j)/5.0
2949 print *,counter,molecule
2950 if (counter.eq.5) then
2952 nres_molec(molecule)=nres_molec(molecule)+1
2955 ! sccor(j,iii)=c(j,ires)
2959 ! print *, "ATOM",atom(1:3)
2960 else if (atom.eq."C5'") then
2961 read (card(19:19),'(a1)') sugar
2962 isugar=sugarcode(sugar,ires)
2967 ! print *,ires,istype(ires)
2970 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2973 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2975 ! write (*,*) card(23:27),ires,itype(ires,1)
2976 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2977 atom.ne.'N' .and. atom.ne.'C' .and. &
2978 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2979 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2980 .and. atom.ne.'P '.and. &
2981 atom(1:1).ne.'H' .and. &
2982 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2984 ! write (iout,*) "sidechain ",atom
2985 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2986 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2987 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2989 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2992 else if ((ions).and.(card(1:6).eq.'HETATM')) then
2993 if (firstion.eq.0) then
2997 dc(j,ires)=sccor(j,iii)
3000 call sccenter(ires,iii,sccor)
3003 read (card(12:16),*) atom
3004 print *,"HETATOM", atom
3005 read (card(18:20),'(a3)') res
3006 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3007 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3008 .or.(atom(1:2).eq.'K ')) &
3011 if (molecule.ne.5) molecprev=molecule
3013 nres_molec(molecule)=nres_molec(molecule)+1
3014 print *,"HERE",nres_molec(molecule)
3015 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
3016 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3020 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3021 if (ires.eq.0) return
3022 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3025 if ((ires_old.ne.ires).and.(molecule.ne.5)) &
3026 nres_molec(molecule)=nres_molec(molecule)-2
3027 ! print *,'I have', nres_molec(:)
3029 do k=1,4 ! ions are without dummy
3030 if (nres_molec(k).eq.0) cycle
3032 ! write (iout,*) i,itype(i,1)
3033 ! if (itype(i,1).eq.ntyp1) then
3034 ! write (iout,*) "dummy",i,itype(i,1)
3036 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3037 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3041 if (itype(i,k).eq.ntyp1_molec(k)) then
3042 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3043 if (itype(i+2,k).eq.0) then
3044 ! print *,"masz sieczke"
3046 if (itype(i+2,j).ne.0) then
3048 itype(i+1,j)=ntyp1_molec(j)
3049 nres_molec(k)=nres_molec(k)-1
3050 nres_molec(j)=nres_molec(j)+1
3056 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3057 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3058 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3059 ! if (unres_pdb) then
3060 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3061 ! print *,i,'tu dochodze'
3062 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3070 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3074 dcj=(c(j,i-2)-c(j,i-3))/2.0
3075 if (dcj.eq.0) dcj=1.23591524223
3080 else !itype(i+1,1).eq.ntyp1
3081 ! if (unres_pdb) then
3082 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3083 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3090 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3094 dcj=(c(j,i+3)-c(j,i+2))/2.0
3095 if (dcj.eq.0) dcj=1.23591524223
3100 endif !itype(i+1,1).eq.ntyp1
3101 endif !itype.eq.ntyp1
3105 ! Calculate the CM of the last side chain.
3109 dc(j,ires)=sccor(j,iii)
3112 call sccenter(ires,iii,sccor)
3118 ! print *,"molecule",molecule
3119 if ((itype(nres,1).ne.10)) then
3121 if (molecule.eq.5) molecule=molecprev
3122 itype(nres,molecule)=ntyp1_molec(molecule)
3123 nres_molec(molecule)=nres_molec(molecule)+1
3124 ! if (unres_pdb) then
3125 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3126 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3133 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3137 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3138 c(j,nres)=c(j,nres-1)+dcj
3139 c(j,2*nres)=c(j,nres)
3143 ! print *,'I have',nres, nres_molec(:)
3145 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3147 if (nres.ne.nres0) then
3148 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3150 stop "Error nres value in WHAM input"
3153 !---------------------------------
3154 !el reallocate tables
3157 ! hfrag_alloc(j,i)=hfrag(j,i)
3160 ! bfrag_alloc(j,i)=bfrag(j,i)
3166 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3167 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3168 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3169 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3173 ! hfrag(j,i)=hfrag_alloc(j,i)
3178 ! bfrag(j,i)=bfrag_alloc(j,i)
3181 !el end reallocate tables
3182 !---------------------------------
3190 c(j,2*nres)=c(j,nres)
3193 if (itype(1,1).eq.ntyp1) then
3197 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3198 call refsys(2,3,4,e1,e2,e3,fail)
3205 c(j,1)=c(j,2)-1.9d0*e2(j)
3209 dcj=(c(j,4)-c(j,3))/2.0
3215 ! First lets assign correct dummy to correct type of chain
3217 if (itype(1,1).eq.ntyp1) then
3218 if (itype(2,1).eq.0) then
3220 if (itype(2,j).ne.0) then
3222 itype(1,j)=ntyp1_molec(j)
3223 nres_molec(1)=nres_molec(1)-1
3224 nres_molec(j)=nres_molec(j)+1
3231 print *,'I have',nres, nres_molec(:)
3233 ! Copy the coordinates to reference coordinates
3239 ! Calculate internal coordinates.
3241 write (iout,'(/a)') &
3242 "Cartesian coordinates of the reference structure"
3243 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3244 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3246 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3247 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3248 (c(j,ires+nres),j=1,3)
3251 ! znamy już nres wiec mozna alokowac tablice
3252 ! Calculate internal coordinates.
3253 if(me.eq.king.or..not.out1file)then
3254 write (iout,'(a)') &
3255 "Backbone and SC coordinates as read from the PDB"
3257 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3258 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3259 (c(j,nres+ires),j=1,3)
3262 ! NOW LETS ROCK! SORTING
3263 allocate(c_temporary(3,2*nres))
3264 allocate(itype_temporary(nres,5))
3265 allocate(molnum(nres+1))
3266 allocate(istype_temp(nres))
3267 itype_temporary(:,:)=0
3271 if (itype(i,k).ne.0) then
3273 c_temporary(j,seqalingbegin)=c(j,i)
3274 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3277 itype_temporary(seqalingbegin,k)=itype(i,k)
3278 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3279 istype_temp(seqalingbegin)=istype(i)
3280 molnum(seqalingbegin)=k
3281 seqalingbegin=seqalingbegin+1
3287 c(j,i)=c_temporary(j,i)
3292 itype(i,k)=itype_temporary(i,k)
3293 istype(i)=istype_temp(i)
3296 if (itype(1,1).eq.ntyp1) then
3300 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3301 call refsys(2,3,4,e1,e2,e3,fail)
3308 c(j,1)=c(j,2)-1.9d0*e2(j)
3312 dcj=(c(j,4)-c(j,3))/2.0
3320 write (iout,'(/a)') &
3321 "Cartesian coordinates of the reference structure after sorting"
3322 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3323 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3325 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3326 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3327 (c(j,ires+nres),j=1,3)
3331 ! print *,seqalingbegin,nres
3332 if(.not.allocated(vbld)) then
3333 allocate(vbld(2*nres))
3338 if(.not.allocated(vbld_inv)) then
3339 allocate(vbld_inv(2*nres))
3345 if(.not.allocated(theta)) then
3346 allocate(theta(nres+2))
3350 if(.not.allocated(phi)) allocate(phi(nres+2))
3351 if(.not.allocated(alph)) allocate(alph(nres+2))
3352 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3353 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3354 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3355 if(.not.allocated(costtab)) allocate(costtab(nres))
3356 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3357 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3358 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3359 if(.not.allocated(xxref)) allocate(xxref(nres))
3360 if(.not.allocated(yyref)) allocate(yyref(nres))
3361 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3362 if(.not.allocated(dc_norm)) then
3363 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3364 allocate(dc_norm(3,0:2*nres+2))
3368 call int_from_cart(.true.,.false.)
3369 call sc_loc_geom(.false.)
3371 thetaref(i)=theta(i)
3381 dc(j,i)=c(j,i+1)-c(j,i)
3382 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3387 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3388 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3390 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3394 ! Copy the coordinates to reference coordinates
3395 ! Splits to single chain if occurs
3399 ! cref(j,i,cou)=c(j,i)
3403 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3404 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3405 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3406 !-----------------------------
3410 write (iout,*) "symetr", symetr
3413 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3415 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3418 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3423 cref(j,i,cou)=c(j,i)
3424 cref(j,i+nres,cou)=c(j,i+nres)
3426 chain_rep(j,lll,kkk)=c(j,i)
3427 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3431 write (iout,*) chain_length
3432 if (chain_length.eq.0) chain_length=nres
3434 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3435 chain_rep(j,chain_length+nres,symetr) &
3436 =chain_rep(j,chain_length+nres,1)
3439 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3441 ! do kkk=1,chain_length
3442 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3446 ! makes copy of chains
3447 write (iout,*) "symetr", symetr
3452 if (symetr.gt.1) then
3459 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3465 ! write (iout,*) i,icha
3466 do lll=1,chain_length
3468 if (cou.le.nres) then
3470 kupa=mod(lll,chain_length)
3471 iprzes=(kkk-1)*chain_length+lll
3472 if (kupa.eq.0) kupa=chain_length
3473 ! write (iout,*) "kupa", kupa
3474 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3475 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3482 !-koniec robienia kopii
3485 write (iout,*) "nowa struktura", nperm
3487 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3489 cref(3,i,kkk),cref(1,nres+i,kkk),&
3490 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3492 100 format (//' alpha-carbon coordinates ',&
3493 ' centroid coordinates'/ &
3494 ' ', 6X,'X',11X,'Y',11X,'Z', &
3495 10X,'X',11X,'Y',11X,'Z')
3496 110 format (a,'(',i3,')',6f12.5)
3502 bfrag(i,j)=bfrag(i,j)-ishift
3508 hfrag(i,j)=hfrag(i,j)-ishift
3514 end subroutine readpdb
3515 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3516 !-----------------------------------------------------------------------------
3518 !-----------------------------------------------------------------------------
3519 subroutine read_control
3533 use random, only: random_init
3534 ! implicit real*8 (a-h,o-z)
3535 ! include 'DIMENSIONS'
3537 use prng, only:prng_restart
3539 logical :: OKRandom!, prng_restart
3542 ! include 'COMMON.IOUNITS'
3543 ! include 'COMMON.TIME1'
3544 ! include 'COMMON.THREAD'
3545 ! include 'COMMON.SBRIDGE'
3546 ! include 'COMMON.CONTROL'
3547 ! include 'COMMON.MCM'
3548 ! include 'COMMON.MAP'
3549 ! include 'COMMON.HEADER'
3550 ! include 'COMMON.CSA'
3551 ! include 'COMMON.CHAIN'
3552 ! include 'COMMON.MUCA'
3553 ! include 'COMMON.MD'
3554 ! include 'COMMON.FFIELD'
3555 ! include 'COMMON.INTERACT'
3556 ! include 'COMMON.SETUP'
3557 !el integer :: KDIAG,ICORFL,IXDR
3558 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3559 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3560 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3561 ! character(len=80) :: ucase
3562 character(len=640) :: controlcard
3564 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3570 read (INP,'(a)') titel
3571 call card_concat(controlcard,.true.)
3572 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3573 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3574 call reada(controlcard,'SEED',seed,0.0D0)
3575 call random_init(seed)
3576 ! Set up the time limit (caution! The time must be input in minutes!)
3577 read_cart=index(controlcard,'READ_CART').gt.0
3578 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3579 call readi(controlcard,'SYM',symetr,1)
3580 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3581 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3582 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3583 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3584 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3585 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3586 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3587 call reada(controlcard,'DRMS',drms,0.1D0)
3588 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3589 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3590 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3591 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3592 write (iout,'(a,f10.1)')'DRMS = ',drms
3593 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3594 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3596 call readi(controlcard,'NZ_START',nz_start,0)
3597 call readi(controlcard,'NZ_END',nz_end,0)
3598 call readi(controlcard,'IZ_SC',iz_sc,0)
3599 timlim=60.0D0*timlim
3600 safety = 60.0d0*safety
3603 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3604 !C SHIELD keyword sets if the shielding effect of side-chains is used
3605 !C 0 denots no shielding is used all peptide are equally despite the
3606 !C solvent accesible area
3607 !C 1 the newly introduced function
3608 !C 2 reseved for further possible developement
3609 call readi(controlcard,'SHIELD',shield_mode,0)
3610 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3611 write(iout,*) "shield_mode",shield_mode
3612 !C Varibles set size of box
3613 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3614 protein=index(controlcard,"PROTEIN").gt.0
3615 ions=index(controlcard,"IONS").gt.0
3616 nucleic=index(controlcard,"NUCLEIC").gt.0
3617 write (iout,*) "with_theta_constr ",with_theta_constr
3618 AFMlog=(index(controlcard,'AFM'))
3619 selfguide=(index(controlcard,'SELFGUIDE'))
3620 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3621 call readi(controlcard,'GENCONSTR',genconstr,0)
3622 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3623 call reada(controlcard,'BOXY',boxysize,100.0d0)
3624 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3625 call readi(controlcard,'TUBEMOD',tubemode,0)
3626 write (iout,*) TUBEmode,"TUBEMODE"
3627 if (TUBEmode.gt.0) then
3628 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3629 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3630 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3631 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3632 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3633 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3634 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3635 buftubebot=bordtubebot+tubebufthick
3636 buftubetop=bordtubetop-tubebufthick
3639 ! CUTOFFF ON ELECTROSTATICS
3640 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3641 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3642 write(iout,*) "R_CUT_ELE=",r_cut_ele
3643 ! Lipidic parameters
3644 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3645 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3646 if (lipthick.gt.0.0d0) then
3647 bordliptop=(boxzsize+lipthick)/2.0
3648 bordlipbot=bordliptop-lipthick
3649 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3650 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3651 buflipbot=bordlipbot+lipbufthick
3652 bufliptop=bordliptop-lipbufthick
3653 if ((lipbufthick*2.0d0).gt.lipthick) &
3654 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3655 endif !lipthick.gt.0
3656 write(iout,*) "bordliptop=",bordliptop
3657 write(iout,*) "bordlipbot=",bordlipbot
3658 write(iout,*) "bufliptop=",bufliptop
3659 write(iout,*) "buflipbot=",buflipbot
3660 write (iout,*) "SHIELD MODE",shield_mode
3662 !C-------------------------
3663 minim=(index(controlcard,'MINIMIZE').gt.0)
3664 dccart=(index(controlcard,'CART').gt.0)
3665 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3666 overlapsc=.not.overlapsc
3667 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3668 searchsc=.not.searchsc
3669 sideadd=(index(controlcard,'SIDEADD').gt.0)
3670 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3671 outpdb=(index(controlcard,'PDBOUT').gt.0)
3672 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3673 pdbref=(index(controlcard,'PDBREF').gt.0)
3674 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3675 indpdb=index(controlcard,'PDBSTART')
3676 extconf=(index(controlcard,'EXTCONF').gt.0)
3677 call readi(controlcard,'IPRINT',iprint,0)
3678 call readi(controlcard,'MAXGEN',maxgen,10000)
3679 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3680 call readi(controlcard,"KDIAG",kdiag,0)
3681 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3682 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3683 write (iout,*) "RESCALE_MODE",rescale_mode
3684 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3685 if (index(controlcard,'REGULAR').gt.0.0D0) then
3686 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3690 if (index(controlcard,'CHECKGRAD').gt.0) then
3692 if (index(controlcard,'CART').gt.0) then
3694 elseif (index(controlcard,'CARINT').gt.0) then
3699 elseif (index(controlcard,'THREAD').gt.0) then
3701 call readi(controlcard,'THREAD',nthread,0)
3702 if (nthread.gt.0) then
3703 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3706 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3707 stop 'Error termination in Read_Control.'
3709 else if (index(controlcard,'MCMA').gt.0) then
3711 else if (index(controlcard,'MCEE').gt.0) then
3713 else if (index(controlcard,'MULTCONF').gt.0) then
3715 else if (index(controlcard,'MAP').gt.0) then
3717 call readi(controlcard,'MAP',nmap,0)
3718 else if (index(controlcard,'CSA').gt.0) then
3720 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3722 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3725 !fcm else if (index(controlcard,'MCMF').gt.0) then
3727 else if (index(controlcard,'SOFTREG').gt.0) then
3729 else if (index(controlcard,'CHECK_BOND').gt.0) then
3731 else if (index(controlcard,'TEST').gt.0) then
3733 else if (index(controlcard,'MD').gt.0) then
3735 else if (index(controlcard,'RE ').gt.0) then
3739 lmuca=index(controlcard,'MUCA').gt.0
3740 call readi(controlcard,'MUCADYN',mucadyn,0)
3741 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3742 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3744 write (iout,*) 'MUCADYN=',mucadyn
3745 write (iout,*) 'MUCASMOOTH=',muca_smooth
3748 iscode=index(controlcard,'ONE_LETTER')
3749 indphi=index(controlcard,'PHI')
3750 indback=index(controlcard,'BACK')
3751 iranconf=index(controlcard,'RAND_CONF')
3752 i2ndstr=index(controlcard,'USE_SEC_PRED')
3753 gradout=index(controlcard,'GRADOUT').gt.0
3754 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3755 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3756 if (me.eq.king .or. .not.out1file ) &
3757 write (iout,*) "DISTCHAINMAX",distchainmax
3759 if(me.eq.king.or..not.out1file) &
3760 write (iout,'(2a)') diagmeth(kdiag),&
3761 ' routine used to diagonalize matrices.'
3762 if (shield_mode.gt.0) then
3764 !C VSolvSphere the volume of solving sphere
3766 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3767 !C there will be no distinction between proline peptide group and normal peptide
3768 !C group in case of shielding parameters
3769 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3770 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3771 write (iout,*) VSolvSphere,VSolvSphere_div
3772 !C long axis of side chain
3774 long_r_sidechain(i)=vbldsc0(1,i)
3775 short_r_sidechain(i)=sigma0(i)
3776 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3782 end subroutine read_control
3783 !-----------------------------------------------------------------------------
3784 subroutine read_REMDpar
3786 ! Read REMD settings
3793 use control_data, only:out1file
3794 ! implicit real*8 (a-h,o-z)
3795 ! include 'DIMENSIONS'
3796 ! include 'COMMON.IOUNITS'
3797 ! include 'COMMON.TIME1'
3798 ! include 'COMMON.MD'
3801 !el include 'COMMON.LANGEVIN'
3803 !el include 'COMMON.LANGEVIN.lang0'
3805 ! include 'COMMON.INTERACT'
3806 ! include 'COMMON.NAMES'
3807 ! include 'COMMON.GEO'
3808 ! include 'COMMON.REMD'
3809 ! include 'COMMON.CONTROL'
3810 ! include 'COMMON.SETUP'
3811 ! character(len=80) :: ucase
3812 character(len=320) :: controlcard
3813 character(len=3200) :: controlcard1
3814 integer :: iremd_m_total
3817 ! real(kind=8) :: var,ene
3819 if(me.eq.king.or..not.out1file) &
3820 write (iout,*) "REMD setup"
3822 call card_concat(controlcard,.true.)
3823 call readi(controlcard,"NREP",nrep,3)
3824 call readi(controlcard,"NSTEX",nstex,1000)
3825 call reada(controlcard,"RETMIN",retmin,10.0d0)
3826 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3827 mremdsync=(index(controlcard,'SYNC').gt.0)
3828 call readi(controlcard,"NSYN",i_sync_step,100)
3829 restart1file=(index(controlcard,'REST1FILE').gt.0)
3830 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3831 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3832 if(max_cache_traj_use.gt.max_cache_traj) &
3833 max_cache_traj_use=max_cache_traj
3834 if(me.eq.king.or..not.out1file) then
3835 !d if (traj1file) then
3836 !rc caching is in testing - NTWX is not ignored
3837 !d write (iout,*) "NTWX value is ignored"
3838 !d write (iout,*) " trajectory is stored to one file by master"
3839 !d write (iout,*) " before exchange at NSTEX intervals"
3841 write (iout,*) "NREP= ",nrep
3842 write (iout,*) "NSTEX= ",nstex
3843 write (iout,*) "SYNC= ",mremdsync
3844 write (iout,*) "NSYN= ",i_sync_step
3845 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3848 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3849 if (index(controlcard,'TLIST').gt.0) then
3851 call card_concat(controlcard1,.true.)
3852 read(controlcard1,*) (remd_t(i),i=1,nrep)
3853 if(me.eq.king.or..not.out1file) &
3854 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3857 if (index(controlcard,'MLIST').gt.0) then
3859 call card_concat(controlcard1,.true.)
3860 read(controlcard1,*) (remd_m(i),i=1,nrep)
3861 if(me.eq.king.or..not.out1file) then
3862 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3865 iremd_m_total=iremd_m_total+remd_m(i)
3867 write (iout,*) 'Total number of replicas ',iremd_m_total
3870 if(me.eq.king.or..not.out1file) &
3871 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3873 end subroutine read_REMDpar
3874 !-----------------------------------------------------------------------------
3875 subroutine read_MDpar
3879 use control_data, only: r_cut,rlamb,out1file
3881 use geometry_data, only: pi
3883 ! implicit real*8 (a-h,o-z)
3884 ! include 'DIMENSIONS'
3885 ! include 'COMMON.IOUNITS'
3886 ! include 'COMMON.TIME1'
3887 ! include 'COMMON.MD'
3890 !el include 'COMMON.LANGEVIN'
3892 !el include 'COMMON.LANGEVIN.lang0'
3894 ! include 'COMMON.INTERACT'
3895 ! include 'COMMON.NAMES'
3896 ! include 'COMMON.GEO'
3897 ! include 'COMMON.SETUP'
3898 ! include 'COMMON.CONTROL'
3899 ! include 'COMMON.SPLITELE'
3900 ! character(len=80) :: ucase
3901 character(len=320) :: controlcard
3906 call card_concat(controlcard,.true.)
3907 call readi(controlcard,"NSTEP",n_timestep,1000000)
3908 call readi(controlcard,"NTWE",ntwe,100)
3909 call readi(controlcard,"NTWX",ntwx,1000)
3910 call reada(controlcard,"DT",d_time,1.0d-1)
3911 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3912 call reada(controlcard,"DAMAX",damax,1.0d1)
3913 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3914 call readi(controlcard,"LANG",lang,0)
3915 RESPA = index(controlcard,"RESPA") .gt. 0
3916 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3917 ntime_split0=ntime_split
3918 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3919 ntime_split0=ntime_split
3920 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3921 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3922 rest = index(controlcard,"REST").gt.0
3923 tbf = index(controlcard,"TBF").gt.0
3924 usampl = index(controlcard,"USAMPL").gt.0
3925 mdpdb = index(controlcard,"MDPDB").gt.0
3926 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3927 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3928 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3929 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3930 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3931 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3932 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3933 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3934 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3935 large = index(controlcard,"LARGE").gt.0
3936 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3937 rattle = index(controlcard,"RATTLE").gt.0
3938 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3944 if(me.eq.king.or..not.out1file) then
3946 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3948 write (iout,'(a)') "The units are:"
3949 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3950 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3951 " acceleration: angstrom/(48.9 fs)**2"
3952 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3954 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3955 write (iout,'(a60,f10.5,a)') &
3956 "Initial time step of numerical integration:",d_time,&
3958 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3960 write (iout,'(2a,i4,a)') &
3961 "A-MTS algorithm used; initial time step for fast-varying",&
3962 " short-range forces split into",ntime_split," steps."
3963 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3964 r_cut," lambda",rlamb
3966 write (iout,'(2a,f10.5)') &
3967 "Maximum acceleration threshold to reduce the time step",&
3968 "/increase split number:",damax
3969 write (iout,'(2a,f10.5)') &
3970 "Maximum predicted energy drift to reduce the timestep",&
3971 "/increase split number:",edriftmax
3972 write (iout,'(a60,f10.5)') &
3973 "Maximum velocity threshold to reduce velocities:",dvmax
3974 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3975 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3976 if (rattle) write (iout,'(a60)') &
3977 "Rattle algorithm used to constrain the virtual bonds"
3981 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3982 call reada(controlcard,"RWAT",rwat,1.4d0)
3983 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3984 surfarea=index(controlcard,"SURFAREA").gt.0
3985 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3986 if(me.eq.king.or..not.out1file)then
3987 write (iout,'(/a,$)') "Langevin dynamics calculation"
3989 write (iout,'(a/)') &
3990 " with direct integration of Langevin equations"
3991 else if (lang.eq.2) then
3992 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3993 else if (lang.eq.3) then
3994 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3995 else if (lang.eq.4) then
3996 write (iout,'(a/)') " in overdamped mode"
3998 write (iout,'(//a,i5)') &
3999 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4002 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4003 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4004 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4005 write (iout,'(a60,f10.5)') &
4006 "Scaling factor of the friction forces:",scal_fric
4007 if (surfarea) write (iout,'(2a,i10,a)') &
4008 "Friction coefficients will be scaled by solvent-accessible",&
4009 " surface area every",reset_fricmat," steps."
4011 ! Calculate friction coefficients and bounds of stochastic forces
4012 eta=6*pi*cPoise*etawat
4013 if(me.eq.king.or..not.out1file) &
4014 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4017 do j=1,5 !types of molecules
4018 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4019 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4021 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4022 do j=1,5 !types of molecules
4024 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4025 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4029 if(me.eq.king.or..not.out1file)then
4030 write (iout,'(/2a/)') &
4031 "Radii of site types and friction coefficients and std's of",&
4032 " stochastic forces of fully exposed sites"
4033 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4035 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4036 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4040 if(me.eq.king.or..not.out1file)then
4041 write (iout,'(a)') "Berendsen bath calculation"
4042 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4043 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4045 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4046 count_reset_moment," steps"
4048 write (iout,'(a,i10,a)') &
4049 "Velocities will be reset at random every",count_reset_vel,&
4053 if(me.eq.king.or..not.out1file) &
4054 write (iout,'(a31)') "Microcanonical mode calculation"
4056 if(me.eq.king.or..not.out1file)then
4057 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4059 write(iout,*) "MD running with constraints."
4060 write(iout,*) "Equilibration time ", eq_time, " mtus."
4061 write(iout,*) "Constraining ", nfrag," fragments."
4062 write(iout,*) "Length of each fragment, weight and q0:"
4064 write (iout,*) "Set of restraints #",iset
4066 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4067 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4069 write(iout,*) "constraints between ", npair, "fragments."
4070 write(iout,*) "constraint pairs, weights and q0:"
4072 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4073 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4075 write(iout,*) "angle constraints within ", nfrag_back,&
4076 "backbone fragments."
4077 write(iout,*) "fragment, weights:"
4079 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4080 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4081 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4084 iset=mod(kolor,nset)+1
4087 if(me.eq.king.or..not.out1file) &
4088 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4090 end subroutine read_MDpar
4091 !-----------------------------------------------------------------------------
4095 ! implicit real*8 (a-h,o-z)
4096 ! include 'DIMENSIONS'
4097 ! include 'COMMON.MAP'
4098 ! include 'COMMON.IOUNITS'
4099 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4100 character(len=80) :: mapcard !,ucase
4103 ! real(kind=8) :: var,ene
4106 read (inp,'(a)') mapcard
4107 mapcard=ucase(mapcard)
4108 if (index(mapcard,'PHI').gt.0) then
4110 else if (index(mapcard,'THE').gt.0) then
4112 else if (index(mapcard,'ALP').gt.0) then
4114 else if (index(mapcard,'OME').gt.0) then
4117 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4118 stop 'Error - illegal variable spec in MAP card.'
4120 call readi (mapcard,'RES1',res1(imap),0)
4121 call readi (mapcard,'RES2',res2(imap),0)
4122 if (res1(imap).eq.0) then
4123 res1(imap)=res2(imap)
4124 else if (res2(imap).eq.0) then
4125 res2(imap)=res1(imap)
4127 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4128 write (iout,'(a)') &
4129 'Error - illegal definition of variable group in MAP.'
4130 stop 'Error - illegal definition of variable group in MAP.'
4132 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4133 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4134 call readi(mapcard,'NSTEP',nstep(imap),0)
4135 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4136 write (iout,'(a)') &
4137 'Illegal boundary and/or step size specification in MAP.'
4138 stop 'Illegal boundary and/or step size specification in MAP.'
4142 end subroutine map_read
4143 !-----------------------------------------------------------------------------
4146 use control_data, only: vdisulf
4148 ! implicit real*8 (a-h,o-z)
4149 ! include 'DIMENSIONS'
4150 ! include 'COMMON.IOUNITS'
4151 ! include 'COMMON.GEO'
4152 ! include 'COMMON.CSA'
4153 ! include 'COMMON.BANK'
4154 ! include 'COMMON.CONTROL'
4155 ! character(len=80) :: ucase
4156 character(len=620) :: mcmcard
4158 ! integer :: ntf,ik,iw_pdb
4159 ! real(kind=8) :: var,ene
4161 call card_concat(mcmcard,.true.)
4163 call readi(mcmcard,'NCONF',nconf,50)
4164 call readi(mcmcard,'NADD',nadd,0)
4165 call readi(mcmcard,'JSTART',jstart,1)
4166 call readi(mcmcard,'JEND',jend,1)
4167 call readi(mcmcard,'NSTMAX',nstmax,500000)
4168 call readi(mcmcard,'N0',n0,1)
4169 call readi(mcmcard,'N1',n1,6)
4170 call readi(mcmcard,'N2',n2,4)
4171 call readi(mcmcard,'N3',n3,0)
4172 call readi(mcmcard,'N4',n4,0)
4173 call readi(mcmcard,'N5',n5,0)
4174 call readi(mcmcard,'N6',n6,10)
4175 call readi(mcmcard,'N7',n7,0)
4176 call readi(mcmcard,'N8',n8,0)
4177 call readi(mcmcard,'N9',n9,0)
4178 call readi(mcmcard,'N14',n14,0)
4179 call readi(mcmcard,'N15',n15,0)
4180 call readi(mcmcard,'N16',n16,0)
4181 call readi(mcmcard,'N17',n17,0)
4182 call readi(mcmcard,'N18',n18,0)
4184 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4186 call readi(mcmcard,'NDIFF',ndiff,2)
4187 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4188 call readi(mcmcard,'IS1',is1,1)
4189 call readi(mcmcard,'IS2',is2,8)
4190 call readi(mcmcard,'NRAN0',nran0,4)
4191 call readi(mcmcard,'NRAN1',nran1,2)
4192 call readi(mcmcard,'IRR',irr,1)
4193 call readi(mcmcard,'NSEED',nseed,20)
4194 call readi(mcmcard,'NTOTAL',ntotal,10000)
4195 call reada(mcmcard,'CUT1',cut1,2.0d0)
4196 call reada(mcmcard,'CUT2',cut2,5.0d0)
4197 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4198 call readi(mcmcard,'ICMAX',icmax,3)
4199 call readi(mcmcard,'IRESTART',irestart,0)
4200 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4203 call reada(mcmcard,'DELE',dele,20.0d0)
4204 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4205 call readi(mcmcard,'IREF',iref,0)
4206 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4207 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4208 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4209 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4210 write (iout,*) "NCONF_IN",nconf_in
4212 end subroutine csaread
4213 !-----------------------------------------------------------------------------
4217 use control_data, only: MaxMoveType
4220 ! implicit real*8 (a-h,o-z)
4221 ! include 'DIMENSIONS'
4222 ! include 'COMMON.MCM'
4223 ! include 'COMMON.MCE'
4224 ! include 'COMMON.IOUNITS'
4225 ! character(len=80) :: ucase
4226 character(len=320) :: mcmcard
4229 ! real(kind=8) :: var,ene
4231 call card_concat(mcmcard,.true.)
4232 call readi(mcmcard,'MAXACC',maxacc,100)
4233 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4234 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4235 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4236 call readi(mcmcard,'MAXREPM',maxrepm,200)
4237 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4238 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4239 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4240 call reada(mcmcard,'E_UP',e_up,5.0D0)
4241 call reada(mcmcard,'DELTE',delte,0.1D0)
4242 call readi(mcmcard,'NSWEEP',nsweep,5)
4243 call readi(mcmcard,'NSTEPH',nsteph,0)
4244 call readi(mcmcard,'NSTEPC',nstepc,0)
4245 call reada(mcmcard,'TMIN',tmin,298.0D0)
4246 call reada(mcmcard,'TMAX',tmax,298.0D0)
4247 call readi(mcmcard,'NWINDOW',nwindow,0)
4248 call readi(mcmcard,'PRINT_MC',print_mc,0)
4249 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4250 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4251 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4252 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4253 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4254 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4255 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4256 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4257 if (nwindow.gt.0) then
4258 allocate(winstart(nwindow)) !!el (maxres)
4259 allocate(winend(nwindow)) !!el
4260 allocate(winlen(nwindow)) !!el
4261 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4263 winlen(i)=winend(i)-winstart(i)+1
4266 if (tmax.lt.tmin) tmax=tmin
4267 if (tmax.eq.tmin) then
4271 if (nstepc.gt.0 .and. nsteph.gt.0) then
4272 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4273 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4275 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4276 ! Probabilities of different move types
4277 sumpro_type(0)=0.0D0
4278 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4279 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4280 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4281 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4282 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4283 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4284 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4286 print *,'i',i,' sumprotype',sumpro_type(i)
4287 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4288 print *,'i',i,' sumprotype',sumpro_type(i)
4291 end subroutine mcmread
4292 !-----------------------------------------------------------------------------
4293 subroutine read_minim
4296 ! implicit real*8 (a-h,o-z)
4297 ! include 'DIMENSIONS'
4298 ! include 'COMMON.MINIM'
4299 ! include 'COMMON.IOUNITS'
4300 ! character(len=80) :: ucase
4301 character(len=320) :: minimcard
4303 ! integer :: ntf,ik,iw_pdb
4304 ! real(kind=8) :: var,ene
4306 call card_concat(minimcard,.true.)
4307 call readi(minimcard,'MAXMIN',maxmin,2000)
4308 call readi(minimcard,'MAXFUN',maxfun,5000)
4309 call readi(minimcard,'MINMIN',minmin,maxmin)
4310 call readi(minimcard,'MINFUN',minfun,maxmin)
4311 call reada(minimcard,'TOLF',tolf,1.0D-2)
4312 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4313 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4314 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4315 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4316 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4317 'Options in energy minimization:'
4318 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4319 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4320 'MinMin:',MinMin,' MinFun:',MinFun,&
4321 ' TolF:',TolF,' RTolF:',RTolF
4323 end subroutine read_minim
4324 !-----------------------------------------------------------------------------
4325 subroutine openunits
4327 use MD_data, only: usampl
4330 use control_data, only:out1file
4331 use control, only: getenv_loc
4332 ! implicit real*8 (a-h,o-z)
4333 ! include 'DIMENSIONS'
4336 character(len=16) :: form,nodename
4337 integer :: nodelen,ierror,npos
4339 ! include 'COMMON.SETUP'
4340 ! include 'COMMON.IOUNITS'
4341 ! include 'COMMON.MD'
4342 ! include 'COMMON.CONTROL'
4343 integer :: lenpre,lenpot,lentmp !,ilen
4345 character(len=3) :: out1file_text !,ucase
4346 character(len=3) :: ll
4349 ! integer :: ntf,ik,iw_pdb
4350 ! real(kind=8) :: var,ene
4352 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4353 call getenv_loc("PREFIX",prefix)
4355 call getenv_loc("POT",pot)
4356 call getenv_loc("DIRTMP",tmpdir)
4357 call getenv_loc("CURDIR",curdir)
4358 call getenv_loc("OUT1FILE",out1file_text)
4359 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4360 out1file_text=ucase(out1file_text)
4361 if (out1file_text(1:1).eq."Y") then
4364 out1file=fg_rank.gt.0
4369 if (lentmp.gt.0) then
4370 write (*,'(80(1h!))')
4371 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4372 write (*,'(80(1h!))')
4373 write (*,*)"All output files will be on node /tmp directory."
4375 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4376 if (me.eq.king) then
4377 write (*,*) "The master node is ",nodename
4378 else if (fg_rank.eq.0) then
4379 write (*,*) "I am the CG slave node ",nodename
4381 write (*,*) "I am the FG slave node ",nodename
4384 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4385 lenpre = lentmp+lenpre+1
4387 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4388 ! Get the names and open the input files
4389 #if defined(WINIFL) || defined(WINPGI)
4390 open(1,file=pref_orig(:ilen(pref_orig))// &
4391 '.inp',status='old',readonly,shared)
4392 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4393 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4394 ! Get parameter filenames and open the parameter files.
4395 call getenv_loc('BONDPAR',bondname)
4396 open (ibond,file=bondname,status='old',readonly,shared)
4397 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4398 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4399 call getenv_loc('THETPAR',thetname)
4400 open (ithep,file=thetname,status='old',readonly,shared)
4401 call getenv_loc('ROTPAR',rotname)
4402 open (irotam,file=rotname,status='old',readonly,shared)
4403 call getenv_loc('TORPAR',torname)
4404 open (itorp,file=torname,status='old',readonly,shared)
4405 call getenv_loc('TORDPAR',tordname)
4406 open (itordp,file=tordname,status='old',readonly,shared)
4407 call getenv_loc('FOURIER',fouriername)
4408 open (ifourier,file=fouriername,status='old',readonly,shared)
4409 call getenv_loc('ELEPAR',elename)
4410 open (ielep,file=elename,status='old',readonly,shared)
4411 call getenv_loc('SIDEPAR',sidename)
4412 open (isidep,file=sidename,status='old',readonly,shared)
4414 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4415 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4416 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4417 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4418 call getenv_loc('TORPAR_NUCL',torname_nucl)
4419 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4420 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4421 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4422 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4423 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4426 #elif (defined CRAY) || (defined AIX)
4427 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4429 ! print *,"Processor",myrank," opened file 1"
4430 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4431 ! print *,"Processor",myrank," opened file 9"
4432 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4433 ! Get parameter filenames and open the parameter files.
4434 call getenv_loc('BONDPAR',bondname)
4435 open (ibond,file=bondname,status='old',action='read')
4436 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4437 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4439 ! print *,"Processor",myrank," opened file IBOND"
4440 call getenv_loc('THETPAR',thetname)
4441 open (ithep,file=thetname,status='old',action='read')
4442 ! print *,"Processor",myrank," opened file ITHEP"
4443 call getenv_loc('ROTPAR',rotname)
4444 open (irotam,file=rotname,status='old',action='read')
4445 ! print *,"Processor",myrank," opened file IROTAM"
4446 call getenv_loc('TORPAR',torname)
4447 open (itorp,file=torname,status='old',action='read')
4448 ! print *,"Processor",myrank," opened file ITORP"
4449 call getenv_loc('TORDPAR',tordname)
4450 open (itordp,file=tordname,status='old',action='read')
4451 ! print *,"Processor",myrank," opened file ITORDP"
4452 call getenv_loc('SCCORPAR',sccorname)
4453 open (isccor,file=sccorname,status='old',action='read')
4454 ! print *,"Processor",myrank," opened file ISCCOR"
4455 call getenv_loc('FOURIER',fouriername)
4456 open (ifourier,file=fouriername,status='old',action='read')
4457 ! print *,"Processor",myrank," opened file IFOURIER"
4458 call getenv_loc('ELEPAR',elename)
4459 open (ielep,file=elename,status='old',action='read')
4460 ! print *,"Processor",myrank," opened file IELEP"
4461 call getenv_loc('SIDEPAR',sidename)
4462 open (isidep,file=sidename,status='old',action='read')
4464 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4465 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4466 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4467 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4468 call getenv_loc('TORPAR_NUCL',torname_nucl)
4469 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4470 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4471 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4472 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4473 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4475 call getenv_loc('LIPTRANPAR',liptranname)
4476 open (iliptranpar,file=liptranname,status='old',action='read')
4477 call getenv_loc('TUBEPAR',tubename)
4478 open (itube,file=tubename,status='old',action='read')
4479 call getenv_loc('IONPAR',ionname)
4480 open (iion,file=ionname,status='old',action='read')
4482 ! print *,"Processor",myrank," opened file ISIDEP"
4483 ! print *,"Processor",myrank," opened parameter files"
4485 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4486 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4487 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4488 ! Get parameter filenames and open the parameter files.
4489 call getenv_loc('BONDPAR',bondname)
4490 open (ibond,file=bondname,status='old')
4491 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4492 open (ibond_nucl,file=bondname_nucl,status='old')
4494 call getenv_loc('THETPAR',thetname)
4495 open (ithep,file=thetname,status='old')
4496 call getenv_loc('ROTPAR',rotname)
4497 open (irotam,file=rotname,status='old')
4498 call getenv_loc('TORPAR',torname)
4499 open (itorp,file=torname,status='old')
4500 call getenv_loc('TORDPAR',tordname)
4501 open (itordp,file=tordname,status='old')
4502 call getenv_loc('SCCORPAR',sccorname)
4503 open (isccor,file=sccorname,status='old')
4504 call getenv_loc('FOURIER',fouriername)
4505 open (ifourier,file=fouriername,status='old')
4506 call getenv_loc('ELEPAR',elename)
4507 open (ielep,file=elename,status='old')
4508 call getenv_loc('SIDEPAR',sidename)
4509 open (isidep,file=sidename,status='old')
4511 open (ithep_nucl,file=thetname_nucl,status='old')
4512 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4513 open (irotam_nucl,file=rotname_nucl,status='old')
4514 call getenv_loc('TORPAR_NUCL',torname_nucl)
4515 open (itorp_nucl,file=torname_nucl,status='old')
4516 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4517 open (itordp_nucl,file=tordname_nucl,status='old')
4518 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4519 open (isidep_nucl,file=sidename_nucl,status='old')
4521 call getenv_loc('LIPTRANPAR',liptranname)
4522 open (iliptranpar,file=liptranname,status='old')
4523 call getenv_loc('TUBEPAR',tubename)
4524 open (itube,file=tubename,status='old')
4525 call getenv_loc('IONPAR',ionname)
4526 open (iion,file=ionname,status='old')
4528 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4530 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4531 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4532 ! Get parameter filenames and open the parameter files.
4533 call getenv_loc('BONDPAR',bondname)
4534 open (ibond,file=bondname,status='old',action='read')
4535 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4536 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4537 call getenv_loc('THETPAR',thetname)
4538 open (ithep,file=thetname,status='old',action='read')
4539 call getenv_loc('ROTPAR',rotname)
4540 open (irotam,file=rotname,status='old',action='read')
4541 call getenv_loc('TORPAR',torname)
4542 open (itorp,file=torname,status='old',action='read')
4543 call getenv_loc('TORDPAR',tordname)
4544 open (itordp,file=tordname,status='old',action='read')
4545 call getenv_loc('SCCORPAR',sccorname)
4546 open (isccor,file=sccorname,status='old',action='read')
4548 call getenv_loc('THETPARPDB',thetname_pdb)
4549 print *,"thetname_pdb ",thetname_pdb
4550 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4551 print *,ithep_pdb," opened"
4553 call getenv_loc('FOURIER',fouriername)
4554 open (ifourier,file=fouriername,status='old',readonly)
4555 call getenv_loc('ELEPAR',elename)
4556 open (ielep,file=elename,status='old',readonly)
4557 call getenv_loc('SIDEPAR',sidename)
4558 open (isidep,file=sidename,status='old',readonly)
4560 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4561 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4562 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4563 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4564 call getenv_loc('TORPAR_NUCL',torname_nucl)
4565 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4566 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4567 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4568 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4569 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4570 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4571 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4572 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4573 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4574 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4575 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4576 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4577 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4580 call getenv_loc('LIPTRANPAR',liptranname)
4581 open (iliptranpar,file=liptranname,status='old',action='read')
4582 call getenv_loc('TUBEPAR',tubename)
4583 open (itube,file=tubename,status='old',action='read')
4584 call getenv_loc('IONPAR',ionname)
4585 open (iion,file=ionname,status='old',action='read')
4588 call getenv_loc('ROTPARPDB',rotname_pdb)
4589 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4592 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4593 #if defined(WINIFL) || defined(WINPGI)
4594 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4595 #elif (defined CRAY) || (defined AIX)
4596 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4598 open (iscpp_nucl,file=scpname_nucl,status='old')
4600 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4605 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4606 ! Use -DOLDSCP to use hard-coded constants instead.
4608 call getenv_loc('SCPPAR',scpname)
4609 #if defined(WINIFL) || defined(WINPGI)
4610 open (iscpp,file=scpname,status='old',readonly,shared)
4611 #elif (defined CRAY) || (defined AIX)
4612 open (iscpp,file=scpname,status='old',action='read')
4614 open (iscpp,file=scpname,status='old')
4616 open (iscpp,file=scpname,status='old',action='read')
4619 call getenv_loc('PATTERN',patname)
4620 #if defined(WINIFL) || defined(WINPGI)
4621 open (icbase,file=patname,status='old',readonly,shared)
4622 #elif (defined CRAY) || (defined AIX)
4623 open (icbase,file=patname,status='old',action='read')
4625 open (icbase,file=patname,status='old')
4627 open (icbase,file=patname,status='old',action='read')
4630 ! Open output file only for CG processes
4631 ! print *,"Processor",myrank," fg_rank",fg_rank
4632 if (fg_rank.eq.0) then
4634 if (nodes.eq.1) then
4637 npos = dlog10(dfloat(nodes-1))+1
4639 if (npos.lt.3) npos=3
4640 write (liczba,'(i1)') npos
4641 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4643 write (liczba,form) me
4644 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4645 liczba(:ilen(liczba))
4646 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4648 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4650 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4651 liczba(:ilen(liczba))//'.mol2'
4652 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4653 liczba(:ilen(liczba))//'.stat'
4655 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4656 //liczba(:ilen(liczba))//'.stat')
4657 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4660 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4661 liczba(:ilen(liczba))//'.const'
4666 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4667 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4668 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4669 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4670 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4672 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4674 rest2name=prefix(:ilen(prefix))//'.rst'
4676 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4679 #if defined(AIX) || defined(PGI)
4680 if (me.eq.king .or. .not. out1file) &
4681 open(iout,file=outname,status='unknown')
4683 if (fg_rank.gt.0) then
4684 write (liczba,'(i3.3)') myrank/nfgtasks
4685 write (ll,'(bz,i3.3)') fg_rank
4686 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4691 open(igeom,file=intname,status='unknown',position='append')
4692 open(ipdb,file=pdbname,status='unknown')
4693 open(imol2,file=mol2name,status='unknown')
4694 open(istat,file=statname,status='unknown',position='append')
4696 !1out open(iout,file=outname,status='unknown')
4699 if (me.eq.king .or. .not.out1file) &
4700 open(iout,file=outname,status='unknown')
4702 if (fg_rank.gt.0) then
4703 write (liczba,'(i3.3)') myrank/nfgtasks
4704 write (ll,'(bz,i3.3)') fg_rank
4705 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4710 open(igeom,file=intname,status='unknown',access='append')
4711 open(ipdb,file=pdbname,status='unknown')
4712 open(imol2,file=mol2name,status='unknown')
4713 open(istat,file=statname,status='unknown',access='append')
4715 !1out open(iout,file=outname,status='unknown')
4718 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4719 csa_seed=prefix(:lenpre)//'.CSA.seed'
4720 csa_history=prefix(:lenpre)//'.CSA.history'
4721 csa_bank=prefix(:lenpre)//'.CSA.bank'
4722 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4723 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4724 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4725 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4726 csa_int=prefix(:lenpre)//'.int'
4727 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4728 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4729 csa_in=prefix(:lenpre)//'.CSA.in'
4730 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4733 write (iout,'(80(1h-))')
4734 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4735 write (iout,'(80(1h-))')
4736 write (iout,*) "Input file : ",&
4737 pref_orig(:ilen(pref_orig))//'.inp'
4738 write (iout,*) "Output file : ",&
4739 outname(:ilen(outname))
4741 write (iout,*) "Sidechain potential file : ",&
4742 sidename(:ilen(sidename))
4744 write (iout,*) "SCp potential file : ",&
4745 scpname(:ilen(scpname))
4747 write (iout,*) "Electrostatic potential file : ",&
4748 elename(:ilen(elename))
4749 write (iout,*) "Cumulant coefficient file : ",&
4750 fouriername(:ilen(fouriername))
4751 write (iout,*) "Torsional parameter file : ",&
4752 torname(:ilen(torname))
4753 write (iout,*) "Double torsional parameter file : ",&
4754 tordname(:ilen(tordname))
4755 write (iout,*) "SCCOR parameter file : ",&
4756 sccorname(:ilen(sccorname))
4757 write (iout,*) "Bond & inertia constant file : ",&
4758 bondname(:ilen(bondname))
4759 write (iout,*) "Bending parameter file : ",&
4760 thetname(:ilen(thetname))
4761 write (iout,*) "Rotamer parameter file : ",&
4762 rotname(:ilen(rotname))
4765 write (iout,*) "Thetpdb parameter file : ",&
4766 thetname_pdb(:ilen(thetname_pdb))
4769 write (iout,*) "Threading database : ",&
4770 patname(:ilen(patname))
4772 write (iout,*)" DIRTMP : ",&
4774 write (iout,'(80(1h-))')
4777 end subroutine openunits
4778 !-----------------------------------------------------------------------------
4781 use geometry_data, only: nres,dc
4783 ! implicit real*8 (a-h,o-z)
4784 ! include 'DIMENSIONS'
4785 ! include 'COMMON.CHAIN'
4786 ! include 'COMMON.IOUNITS'
4787 ! include 'COMMON.MD'
4790 ! real(kind=8) :: var,ene
4792 open(irest2,file=rest2name,status='unknown')
4793 read(irest2,*) totT,EK,potE,totE,t_bath
4796 ! AL 4/17/17: Now reading d_t(0,:) too
4798 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4801 ! AL 4/17/17: Now reading d_c(0,:) too
4803 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4806 read (irest2,*) iset
4810 end subroutine readrst
4811 !-----------------------------------------------------------------------------
4812 subroutine read_fragments
4816 use control_data, only:out1file
4819 ! implicit real*8 (a-h,o-z)
4820 ! include 'DIMENSIONS'
4824 ! include 'COMMON.SETUP'
4825 ! include 'COMMON.CHAIN'
4826 ! include 'COMMON.IOUNITS'
4827 ! include 'COMMON.MD'
4828 ! include 'COMMON.CONTROL'
4831 ! real(kind=8) :: var,ene
4833 read(inp,*) nset,nfrag,npair,nfrag_back
4835 !el from module energy
4836 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4837 if(.not.allocated(wfrag_back)) then
4838 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4839 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4841 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4842 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4844 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4845 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4848 if(me.eq.king.or..not.out1file) &
4849 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4850 " nfrag_back",nfrag_back
4852 read(inp,*) mset(iset)
4854 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4856 if(me.eq.king.or..not.out1file) &
4857 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4858 ifrag(2,i,iset), qinfrag(i,iset)
4861 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4863 if(me.eq.king.or..not.out1file) &
4864 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4865 ipair(2,i,iset), qinpair(i,iset)
4868 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4869 wfrag_back(3,i,iset),&
4870 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4871 if(me.eq.king.or..not.out1file) &
4872 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4873 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4877 end subroutine read_fragments
4878 !-----------------------------------------------------------------------------
4880 !-----------------------------------------------------------------------------
4884 ! implicit real*8 (a-h,o-z)
4885 ! include 'DIMENSIONS'
4886 ! include 'COMMON.CSA'
4887 ! include 'COMMON.BANK'
4888 ! include 'COMMON.IOUNITS'
4890 ! integer :: ntf,ik,iw_pdb
4891 ! real(kind=8) :: var,ene
4893 open(icsa_in,file=csa_in,status="old",err=100)
4894 read(icsa_in,*) nconf
4895 read(icsa_in,*) jstart,jend
4896 read(icsa_in,*) nstmax
4897 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4898 read(icsa_in,*) nran0,nran1,irr
4899 read(icsa_in,*) nseed
4900 read(icsa_in,*) ntotal,cut1,cut2
4901 read(icsa_in,*) estop
4902 read(icsa_in,*) icmax,irestart
4903 read(icsa_in,*) ntbankm,dele,difcut
4904 read(icsa_in,*) iref,rmscut,pnccut
4905 read(icsa_in,*) ndiff
4912 end subroutine csa_read
4913 !-----------------------------------------------------------------------------
4914 subroutine initial_write
4917 ! implicit real*8 (a-h,o-z)
4918 ! include 'DIMENSIONS'
4919 ! include 'COMMON.CSA'
4920 ! include 'COMMON.BANK'
4921 ! include 'COMMON.IOUNITS'
4923 ! integer :: ntf,ik,iw_pdb
4924 ! real(kind=8) :: var,ene
4926 open(icsa_seed,file=csa_seed,status="unknown")
4927 write(icsa_seed,*) "seed"
4929 #if defined(AIX) || defined(PGI)
4930 open(icsa_history,file=csa_history,status="unknown",&
4933 open(icsa_history,file=csa_history,status="unknown",&
4936 write(icsa_history,*) nconf
4937 write(icsa_history,*) jstart,jend
4938 write(icsa_history,*) nstmax
4939 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4940 write(icsa_history,*) nran0,nran1,irr
4941 write(icsa_history,*) nseed
4942 write(icsa_history,*) ntotal,cut1,cut2
4943 write(icsa_history,*) estop
4944 write(icsa_history,*) icmax,irestart
4945 write(icsa_history,*) ntbankm,dele,difcut
4946 write(icsa_history,*) iref,rmscut,pnccut
4947 write(icsa_history,*) ndiff
4949 write(icsa_history,*)
4952 open(icsa_bank1,file=csa_bank1,status="unknown")
4953 write(icsa_bank1,*) 0
4957 end subroutine initial_write
4958 !-----------------------------------------------------------------------------
4959 subroutine restart_write
4962 ! implicit real*8 (a-h,o-z)
4963 ! include 'DIMENSIONS'
4964 ! include 'COMMON.IOUNITS'
4965 ! include 'COMMON.CSA'
4966 ! include 'COMMON.BANK'
4968 ! integer :: ntf,ik,iw_pdb
4969 ! real(kind=8) :: var,ene
4971 #if defined(AIX) || defined(PGI)
4972 open(icsa_history,file=csa_history,position="append")
4974 open(icsa_history,file=csa_history,access="append")
4976 write(icsa_history,*)
4977 write(icsa_history,*) "This is restart"
4978 write(icsa_history,*)
4979 write(icsa_history,*) nconf
4980 write(icsa_history,*) jstart,jend
4981 write(icsa_history,*) nstmax
4982 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4983 write(icsa_history,*) nran0,nran1,irr
4984 write(icsa_history,*) nseed
4985 write(icsa_history,*) ntotal,cut1,cut2
4986 write(icsa_history,*) estop
4987 write(icsa_history,*) icmax,irestart
4988 write(icsa_history,*) ntbankm,dele,difcut
4989 write(icsa_history,*) iref,rmscut,pnccut
4990 write(icsa_history,*) ndiff
4991 write(icsa_history,*)
4992 write(icsa_history,*) "irestart is: ", irestart
4994 write(icsa_history,*)
4998 end subroutine restart_write
4999 !-----------------------------------------------------------------------------
5001 !-----------------------------------------------------------------------------
5002 subroutine write_pdb(npdb,titelloc,ee)
5004 ! implicit real*8 (a-h,o-z)
5005 ! include 'DIMENSIONS'
5006 ! include 'COMMON.IOUNITS'
5007 character(len=50) :: titelloc1
5008 character*(*) :: titelloc
5009 character(len=3) :: zahl
5010 character(len=5) :: liczba5
5012 integer :: npdb !,ilen
5016 ! real(kind=8) :: var,ene
5020 if (npdb.lt.1000) then
5021 call numstr(npdb,zahl)
5022 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5024 if (npdb.lt.10000) then
5025 write(liczba5,'(i1,i4)') 0,npdb
5027 write(liczba5,'(i5)') npdb
5029 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5031 call pdbout(ee,titelloc1,ipdb)
5034 end subroutine write_pdb
5035 !-----------------------------------------------------------------------------
5037 !-----------------------------------------------------------------------------
5038 subroutine write_thread_summary
5039 ! Thread the sequence through a database of known structures
5040 use control_data, only: refstr
5042 use energy_data, only: n_ene_comp
5044 ! implicit real*8 (a-h,o-z)
5045 ! include 'DIMENSIONS'
5047 use MPI_data !include 'COMMON.INFO'
5050 ! include 'COMMON.CONTROL'
5051 ! include 'COMMON.CHAIN'
5052 ! include 'COMMON.DBASE'
5053 ! include 'COMMON.INTERACT'
5054 ! include 'COMMON.VAR'
5055 ! include 'COMMON.THREAD'
5056 ! include 'COMMON.FFIELD'
5057 ! include 'COMMON.SBRIDGE'
5058 ! include 'COMMON.HEADER'
5059 ! include 'COMMON.NAMES'
5060 ! include 'COMMON.IOUNITS'
5061 ! include 'COMMON.TIME1'
5063 integer,dimension(maxthread) :: ip
5064 real(kind=8),dimension(0:n_ene) :: energia
5066 integer :: i,j,ii,jj,ipj,ik,kk,ist
5067 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5069 write (iout,'(30x,a/)') &
5070 ' *********** Summary threading statistics ************'
5071 write (iout,'(a)') 'Initial energies:'
5072 write (iout,'(a4,2x,a12,14a14,3a8)') &
5073 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5074 'RMSnat','NatCONT','NNCONT','RMS'
5075 ! Energy sort patterns
5080 enet=ener(n_ene-1,ip(i))
5083 if (ener(n_ene-1,ip(j)).lt.enet) then
5085 enet=ener(n_ene-1,ip(j))
5097 ist=nres_base(2,ii)+ipatt(2,i)
5099 energia(i)=ener0(kk,i)
5101 etot=ener0(n_ene_comp+1,i)
5102 rmsnat=ener0(n_ene_comp+2,i)
5103 rms=ener0(n_ene_comp+3,i)
5104 frac=ener0(n_ene_comp+4,i)
5105 frac_nn=ener0(n_ene_comp+5,i)
5108 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5109 i,str_nam(ii),ist+1,&
5110 (energia(print_order(kk)),kk=1,nprint_ene),&
5111 etot,rmsnat,frac,frac_nn,rms
5113 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5114 i,str_nam(ii),ist+1,&
5115 (energia(print_order(kk)),kk=1,nprint_ene),etot
5118 write (iout,'(//a)') 'Final energies:'
5119 write (iout,'(a4,2x,a12,17a14,3a8)') &
5120 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5121 'RMSnat','NatCONT','NNCONT','RMS'
5125 ist=nres_base(2,ii)+ipatt(2,i)
5127 energia(kk)=ener(kk,ik)
5129 etot=ener(n_ene_comp+1,i)
5130 rmsnat=ener(n_ene_comp+2,i)
5131 rms=ener(n_ene_comp+3,i)
5132 frac=ener(n_ene_comp+4,i)
5133 frac_nn=ener(n_ene_comp+5,i)
5134 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5135 i,str_nam(ii),ist+1,&
5136 (energia(print_order(kk)),kk=1,nprint_ene),&
5137 etot,rmsnat,frac,frac_nn,rms
5139 write (iout,'(/a/)') 'IEXAM array:'
5140 write (iout,'(i5)') nexcl
5142 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5144 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5145 'Max. time for threading step ',max_time_for_thread,&
5146 'Average time for threading step: ',ave_time_for_thread
5148 end subroutine write_thread_summary
5149 !-----------------------------------------------------------------------------
5150 subroutine write_stat_thread(ithread,ipattern,ist)
5152 use energy_data, only: n_ene_comp
5154 ! implicit real*8 (a-h,o-z)
5155 ! include "DIMENSIONS"
5156 ! include "COMMON.CONTROL"
5157 ! include "COMMON.IOUNITS"
5158 ! include "COMMON.THREAD"
5159 ! include "COMMON.FFIELD"
5160 ! include "COMMON.DBASE"
5161 ! include "COMMON.NAMES"
5162 real(kind=8),dimension(0:n_ene) :: energia
5164 integer :: ithread,ipattern,ist,i
5165 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5167 #if defined(AIX) || defined(PGI)
5168 open(istat,file=statname,position='append')
5170 open(istat,file=statname,access='append')
5173 energia(i)=ener(i,ithread)
5175 etot=ener(n_ene_comp+1,ithread)
5176 rmsnat=ener(n_ene_comp+2,ithread)
5177 rms=ener(n_ene_comp+3,ithread)
5178 frac=ener(n_ene_comp+4,ithread)
5179 frac_nn=ener(n_ene_comp+5,ithread)
5180 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5181 ithread,str_nam(ipattern),ist+1,&
5182 (energia(print_order(i)),i=1,nprint_ene),&
5183 etot,rmsnat,frac,frac_nn,rms
5186 end subroutine write_stat_thread
5187 !-----------------------------------------------------------------------------
5189 !-----------------------------------------------------------------------------
5190 end module io_config